Coverage for pythia/sampler.py: 89%

517 statements  

« prev     ^ index     » next       coverage.py v7.3.1, created at 2023-09-08 17:13 +0000

1""" 

2File: pythia/sampler.py 

3Author: Nando Hegemann 

4Gitlab: https://gitlab.com/Nando-Hegemann 

5Description: Sampler classes for generating in random samples and PDF evaluations. 

6SPDX-License-Identifier: LGPL-3.0-or-later OR Hippocratic-3.0-ECO-MEDIA-MIL 

7""" 

8from typing import Callable 

9from abc import ABC, abstractmethod, abstractproperty 

10import warnings 

11import numpy as np 

12import scipy.stats 

13from scipy.integrate import quad 

14from scipy.special import gamma 

15import pythia as pt 

16 

17 

18class Sampler(ABC): 

19 """Base class for all continuous samplers.""" 

20 

21 # set abstract attributes 

22 domain: np.ndarray 

23 

24 @abstractproperty 

25 def dimension(self) -> int: 

26 """Dimension of the ambient space.""" 

27 raise NotImplementedError 

28 

29 @abstractproperty 

30 def mass(self): 

31 """Mass of the sampler distribution. 

32 

33 The integral of the sampler distribution over the domain of 

34 definition. If the density is normalised this value should be one. 

35 """ 

36 raise NotImplementedError 

37 

38 @abstractproperty 

39 def maximum(self) -> float: 

40 """Maximum of the pdf.""" 

41 raise NotImplementedError 

42 

43 @abstractproperty 

44 def mean(self) -> float | np.ndarray: 

45 """Mean value of the pdf.""" 

46 raise NotImplementedError 

47 

48 @abstractproperty 

49 def cov(self) -> float | np.ndarray: 

50 """(Co)Variance of the pdf.""" 

51 raise NotImplementedError 

52 

53 @abstractmethod 

54 def pdf(self, x: np.ndarray) -> np.ndarray: 

55 """Density of the samplers distribution. 

56 

57 Computes the density of the samplers underlying distribution at the 

58 given points `x`. 

59 

60 Parameters 

61 ---------- 

62 x : array_like of shape (..., D) 

63 list of points or single point. `D` is the objects dimension. 

64 

65 Returns 

66 ------- 

67 : 

68 Density values at the points. 

69 """ 

70 raise NotImplementedError 

71 

72 @abstractmethod 

73 def log_pdf(self, x: np.ndarray) -> np.ndarray: 

74 """Log-density of the samplers distribution. 

75 

76 Computes the log-density of the samplers underlying distribution at the 

77 given points `x`. 

78 

79 Parameters 

80 ---------- 

81 x : array_like of shape (..., D) 

82 list of points or single point. `D` is the objects dimension. 

83 

84 Returns 

85 ------- 

86 : 

87 Log-density values at the points. 

88 """ 

89 raise NotImplementedError 

90 

91 @abstractmethod 

92 def grad_x_log_pdf(self, x: np.ndarray) -> np.ndarray: 

93 """Gradient of log-density of the samplers distribution. 

94 

95 Computes the gradient of the log-density of the samplers underlying 

96 distribution at the given points `x`. 

97 

98 Parameters 

99 ---------- 

100 x : array_like of shape (..., D) 

101 list of points or single point. `D` is the objects dimension. 

102 

103 Returns 

104 ------- 

105 : 

106 Gradient values of the log-density at the points with shape 

107 (..., D). 

108 """ 

109 raise NotImplementedError 

110 

111 @abstractmethod 

112 def hess_x_log_pdf(self, x: np.ndarray) -> np.ndarray: 

113 """Hessian of log-density of the samplers distribution. 

114 

115 Computes the Hessian of the log-density of the samplers underlying 

116 distribution at the given points `x`. 

117 

118 Parameters 

119 ---------- 

120 x : array_like of shape (..., D) 

121 list of points or single point. `D` is the objects dimension. 

122 

123 Returns 

124 ------- 

125 : 

126 Hessian values of the log-density at the points with shape 

127 (..., D, D). 

128 """ 

129 raise NotImplementedError 

130 

131 @abstractmethod 

132 def sample(self, shape: list | tuple | np.ndarray) -> np.ndarray: 

133 """Random values in a given shape. 

134 

135 Create an array of the given shape and populate it with random samples 

136 from the samplers distribution. 

137 

138 Parameters 

139 ---------- 

140 shape : array_like, optional 

141 The dimensions of the returned array, should all be positive. 

142 If no argument is given a single Python float is returned. 

143 

144 Returns 

145 ------- 

146 : 

147 Random values of specified shape. 

148 """ 

149 raise NotImplementedError 

150 

151 

152class UniformSampler(Sampler): 

153 """Sampler for univariate uniformly distributed samples on given domain. 

154 

155 Parameters 

156 ---------- 

157 domain : array_like 

158 Interval of support of distribution. 

159 """ 

160 

161 def __init__(self, domain: list | tuple | np.ndarray) -> None: 

162 """Initiate UniformSampler object.""" 

163 self.domain = np.reshape(np.array(domain), (-1, 2)) 

164 assert self.domain.shape == (1, 2) 

165 self._length = self.domain[0, 1] - self.domain[0, 0] 

166 

167 @property 

168 def dimension(self) -> int: 

169 """Parameter dimension.""" 

170 return self.domain.shape[0] 

171 

172 @property 

173 def mass(self) -> float: 

174 """Mass of the PDF.""" 

175 return 1 

176 

177 @property 

178 def maximum(self) -> float: 

179 """Maximum value of the PDF.""" 

180 return 1 / self._length 

181 

182 @property 

183 def mean(self) -> float: 

184 """Mean value of the distribution.""" 

185 return 0.5 * (self.domain[0, 0] + self.domain[0, 1]) 

186 

187 @property 

188 def cov(self) -> float: 

189 """(Co)Variance of the distribution.""" 

190 return 1 / 12 * (self.domain[0, 1] - self.domain[0, 0]) ** 2 

191 

192 @property 

193 def var(self) -> float: 

194 """Variance of the distribution.""" 

195 return self.cov 

196 

197 @property 

198 def std(self) -> float: 

199 """Standard deviation of the distribution.""" 

200 return np.sqrt(self.var) 

201 

202 def pdf(self, x: np.ndarray) -> np.ndarray: 

203 """Evaluate uniform PDF. 

204 

205 Parameters 

206 ---------- 

207 x : array_like 

208 Evaluation points. 

209 

210 Returns 

211 ------- 

212 : 

213 Values of PDF evaluated in `x`. 

214 """ 

215 return scipy.stats.uniform.pdf(x, loc=self.domain[0, 0], scale=self._length) 

216 

217 def log_pdf(self, x: np.ndarray) -> np.ndarray: 

218 """Evaluate uniform log-PDF. 

219 

220 Parameters 

221 ---------- 

222 x : array_like 

223 Evaluation points. 

224 

225 Returns 

226 ------- 

227 : 

228 Values of log-PDF evaluated in `x`. 

229 """ 

230 return scipy.stats.uniform.logpdf(x, loc=self.domain[0, 0], scale=self._length) 

231 

232 def grad_x_log_pdf(self, x: np.ndarray) -> np.ndarray: 

233 """Evaluate gradient of uniform log-PDF. 

234 

235 Parameters 

236 ---------- 

237 x : array_like 

238 Evaluation points. 

239 

240 Returns 

241 ------- 

242 : 

243 Values of gradient (vector valued) of log-PDF evaluated in `x`. 

244 """ 

245 return np.zeros_like(x) 

246 

247 def hess_x_log_pdf(self, x: np.ndarray) -> np.ndarray: 

248 """Evaluate Hessian of uniform log-PDF. 

249 

250 Parameters 

251 ---------- 

252 x : array_like 

253 Evaluation points. 

254 

255 Returns 

256 ------- 

257 : 

258 Values of Hessian (matrix valued) of log-PDF evaluated in `x`. 

259 """ 

260 return np.zeros_like(x) 

261 

262 def sample(self, shape: int | list | tuple | np.ndarray) -> np.ndarray: 

263 """Draw samples from uniform distribution. 

264 

265 Parameters 

266 ---------- 

267 shape : array_like 

268 Shape of the samples. 

269 

270 Returns 

271 ------- 

272 : 

273 Random samples of specified shape. 

274 """ 

275 return np.random.uniform(self.domain[0, 0], self.domain[0, 1], shape) 

276 

277 

278class NormalSampler(Sampler): 

279 """Sampler for univariate normally distributed samples. 

280 

281 Parameters 

282 ---------- 

283 mean : float 

284 Mean of the Gaussian distribution. 

285 var : float 

286 Variance of the Gaussian distribution. 

287 """ 

288 

289 def __init__(self, mean: float, var: float) -> None: 

290 """Initiate NormalSampler object.""" 

291 self.domain = np.array([-np.inf, np.inf], ndmin=2) # shape is (1,2) 

292 self._mean = mean # work-around to comply with ABC Sampler class 

293 assert var >= 0 

294 self.var = var 

295 

296 @property 

297 def mass(self) -> float: 

298 """Mass of the PDF.""" 

299 return 1 

300 

301 @property 

302 def dimension(self) -> float: 

303 """Dimension of the parameters.""" 

304 return self.domain.shape[0] 

305 

306 @property 

307 def maximum(self) -> float: 

308 """Maximum value of the PDF.""" 

309 return 1 / np.sqrt(2 * np.pi * self.var) 

310 

311 @property 

312 def mean(self) -> float: 

313 """Mean value of the distribution.""" 

314 return self._mean 

315 

316 @property 

317 def cov(self) -> float: 

318 """(Co)Variance of the distribution.""" 

319 return self.var 

320 

321 @property 

322 def std(self) -> float: 

323 """Standard deviation.""" 

324 return np.sqrt(self.var) 

325 

326 def pdf(self, x: np.ndarray) -> np.ndarray: 

327 """Evaluate PDF. 

328 

329 Parameters 

330 ---------- 

331 x : array_like 

332 Evaluation points. 

333 

334 Returns 

335 ------- 

336 : 

337 Values of PDF evaluated in `x`. 

338 """ 

339 return scipy.stats.norm.pdf(x, loc=self.mean, scale=np.sqrt(self.var)) 

340 

341 def log_pdf(self, x: np.ndarray) -> np.ndarray: 

342 """Evaluate log-PDF. 

343 

344 Parameters 

345 ---------- 

346 x : array_like 

347 Evaluation points. 

348 

349 Returns 

350 ------- 

351 : 

352 Values of log-PDF evaluated in `x`. 

353 """ 

354 return scipy.stats.norm.logpdf(x, loc=self.mean, scale=np.sqrt(self.var)) 

355 

356 def grad_x_log_pdf(self, x: np.ndarray) -> np.ndarray: 

357 """Evaluate gradient of log-PDF. 

358 

359 Parameters 

360 ---------- 

361 x : array_like 

362 Evaluation points. 

363 

364 Returns 

365 ------- 

366 : 

367 Values of gradient (vector valued) of log-PDF evaluated in `x`. 

368 """ 

369 return -(x - self.mean) / self.var 

370 

371 def hess_x_log_pdf(self, x: np.ndarray) -> np.ndarray: 

372 """Evaluate Hessian of log-PDF. 

373 

374 Parameters 

375 ---------- 

376 x : array_like 

377 Evaluation points. 

378 

379 Returns 

380 ------- 

381 : 

382 Values of Hessian (matrix valued) of log-PDF evaluated in `x`. 

383 """ 

384 return -1 / self.var * np.ones_like(x) 

385 

386 def sample(self, shape: int | list | tuple | np.ndarray) -> np.ndarray: 

387 """Draw samples from distribution. 

388 

389 Parameters 

390 ---------- 

391 shape : array_like 

392 Shape of the samples. 

393 

394 Returns 

395 ------- 

396 : 

397 Random samples of specified shape. 

398 """ 

399 return np.random.normal(self.mean, np.sqrt(self.var), shape) 

400 

401 

402class GammaSampler(Sampler): 

403 """Sampler for univariate Gamma distributed samples on given domain. 

404 

405 Parameters 

406 ---------- 

407 domain : array_like 

408 Supported domain of distribution. 

409 alpha : float 

410 Parameter for Gamma distribution. 

411 beta : float 

412 Parameter for Gamma distribution. 

413 """ 

414 

415 def __init__( 

416 self, domain: list | tuple | np.ndarray, alpha: float, beta: float 

417 ) -> None: 

418 """Initiate GammaSampler object.""" 

419 self.domain = np.reshape(np.array(domain), (-1, 2)) 

420 assert self.domain.shape == (1, 2) 

421 assert self.domain[0, 1] == np.inf 

422 self.alpha = alpha 

423 self.beta = beta 

424 assert self.alpha > 0 and self.beta > 0 

425 

426 @property 

427 def dimension(self) -> float: 

428 """Dimension of the parameters.""" 

429 return 1 

430 

431 @property 

432 def mass(self) -> float: 

433 """Mass of the PDF.""" 

434 return 1 

435 

436 @property 

437 def maximum(self) -> float: 

438 """Maximum value of the PDF. 

439 

440 The maximum of the Gamma distribution is given by 

441 

442 .. math:: 

443 \\max_{x\\in[a,\\infty)} f(x) = 

444 \\begin{cases} 

445 \\infty & \\mbox{if } 0 < \\alpha < 1\\\\ 

446 \\frac{\\beta^\\alpha}{\\Gamma(\\alpha)} & \\mbox{if } \\alpha = 1\\\\ 

447 \\frac{\\beta^\\alpha}{\\Gamma(\\alpha)} \\Bigl(\\frac{\\alpha-1}{\\beta} \\Bigr)^{\\alpha-1} e^{1-\\alpha} & \\mbox{if } \\alpha > 1\\\\ 

448 \\end{cases} 

449 """ 

450 if self.alpha < 1: 

451 return np.inf 

452 

453 if self.alpha == 1: 

454 return self.beta**self.alpha / gamma(self.alpha) 

455 

456 return ( 

457 self.beta**self.alpha 

458 / gamma(self.alpha) 

459 * ((self.alpha - 1) / self.beta) ** (self.alpha - 1) 

460 * np.exp(-(self.alpha - 1)) 

461 ) 

462 

463 @property 

464 def mean(self) -> float: 

465 """Mean value of the distribution.""" 

466 return self.alpha / self.beta + self.domain[0, 0] 

467 

468 @property 

469 def cov(self) -> float: 

470 """(Co)Variance of the distribution.""" 

471 return self.alpha / self.beta**2 

472 

473 @property 

474 def var(self) -> float: 

475 """Variance of the distribution.""" 

476 return self.cov 

477 

478 @property 

479 def std(self) -> float: 

480 """Standard deviation of the distribution.""" 

481 return np.sqrt(self.var) 

482 

483 def pdf(self, x: np.ndarray) -> np.ndarray: 

484 """Evaluate PDF. 

485 

486 Parameters 

487 ---------- 

488 x : array_like 

489 Evaluation points. 

490 

491 Returns 

492 ------- 

493 : 

494 Values of PDF evaluated in `x`. 

495 """ 

496 y = x - self.domain[0, 0] 

497 return scipy.stats.gamma.pdf(y, a=self.alpha, scale=1.0 / self.beta) 

498 

499 def log_pdf(self, x: np.ndarray) -> np.ndarray: 

500 """Evaluate log-PDF. 

501 

502 Parameters 

503 ---------- 

504 x : array_like 

505 Evaluation points. 

506 

507 Returns 

508 ------- 

509 : 

510 Values of log-PDF evaluated in `x`. 

511 """ 

512 y = x - self.domain[0, 0] 

513 return scipy.stats.gamma.logpdf(y, a=self.alpha, scale=1.0 / self.beta) 

514 

515 def grad_x_log_pdf(self, x: np.ndarray) -> np.ndarray: 

516 """Evaluate gradient of log-PDF. 

517 

518 .. note:: 

519 Not yet implemented. 

520 

521 Parameters 

522 ---------- 

523 x : array_like 

524 Evaluation points. 

525 

526 Returns 

527 ------- 

528 : 

529 Values of gradient (vector valued) of log-PDF evaluated in `x`. 

530 """ 

531 raise NotImplementedError 

532 

533 def hess_x_log_pdf(self, x: np.ndarray) -> np.ndarray: 

534 """Evaluate Hessian of log-PDF. 

535 

536 .. note:: 

537 Not yet implemented. 

538 

539 Parameters 

540 ---------- 

541 x : array_like 

542 Evaluation points. 

543 

544 Returns 

545 ------- 

546 : 

547 Values of Hessian (matrix valued) of log-PDF evaluated in `x`. 

548 """ 

549 raise NotImplementedError 

550 

551 def sample(self, shape: int | list | tuple | np.ndarray) -> np.ndarray: 

552 """Draw samples from distribution. 

553 

554 Parameters 

555 ---------- 

556 shape : array_like 

557 Shape of the samples. 

558 

559 Returns 

560 ------- 

561 : 

562 Random samples of specified shape. 

563 """ 

564 ret = np.random.gamma(self.alpha, 1.0 / self.beta, shape) 

565 return ret + self.domain[0, 0] 

566 

567 

568class BetaSampler(Sampler): 

569 """Sampler for univariate Beta distributed samples on given domain. 

570 

571 Parameters 

572 ---------- 

573 domain : array_like 

574 Supported domain of distribution. 

575 alpha : float 

576 Parameter for Beta distribution. 

577 beta : float 

578 Parameter for Beta distribution. 

579 """ 

580 

581 def __init__( 

582 self, domain: list | tuple | np.ndarray, alpha: float, beta: float 

583 ) -> None: 

584 """Initiate BetaSampler object.""" 

585 self.domain = np.reshape(np.array(domain), (-1, 2)) 

586 assert self.domain.shape == (1, 2) 

587 self.length = self.domain[0, 1] - self.domain[0, 0] 

588 self.alpha = alpha 

589 self.beta = beta 

590 assert self.alpha > 0 and self.beta > 0 

591 

592 @property 

593 def dimension(self): 

594 """Dimension of the parameters.""" 

595 return 1 

596 

597 @property 

598 def mass(self): 

599 """Mass of the PDF.""" 

600 return 1 

601 

602 @property 

603 def maximum(self): 

604 """Maximum value of the PDF. 

605 

606 The maximum of the Beta distribution is given by 

607 

608 .. math:: 

609 \\max_{x\\in[a,b]} f(x) = 

610 \\begin{cases} 

611 \\infty & \\mbox{if } 0 < \\alpha < 1 \\mbox{ or } 0 < \\beta < 1,\\\\ 

612 \\frac{1}{(b-a)B(\\alpha,\\beta)} & \\mbox{if } \\alpha = 1 \\mbox{ or } \\beta = 1,\\\\ 

613 \\frac{(\\alpha-1)^{\\alpha-1}(\\beta-1)^{\\beta-1}}{(\\alpha+\\beta-2)^{\\alpha+\\beta-2}(b-a)B(\\alpha,\\beta)} & \\mbox{if } \\alpha > 1 \\mbox{ and } \\beta > 1, \\\\ 

614 \\end{cases} 

615 

616 where :math:`B(\\alpha,\\beta)` denotes the Beta-function. 

617 """ 

618 if self.alpha < 1 or self.beta < 1: 

619 return np.inf 

620 

621 Beta = gamma(self.alpha) * gamma(self.beta) / gamma(self.alpha + self.beta) 

622 if self.alpha == 1 or self.beta == 1: 

623 return 1 / (Beta * self.length) 

624 

625 val = ( 

626 (self.alpha - 1) ** (self.alpha - 1) 

627 * (self.beta - 1) ** (self.beta - 1) 

628 / (self.alpha + self.beta - 2) ** (self.alpha + self.beta - 2) 

629 ) 

630 return val / (Beta * self.length) 

631 

632 @property 

633 def mean(self) -> float: 

634 """Mean value of the distribution.""" 

635 shift = self.domain[0, 0] 

636 scale = self.length 

637 return scale * self.alpha / (self.alpha + self.beta) + shift 

638 

639 @property 

640 def cov(self) -> float: 

641 """(Co)Variance of the distribution.""" 

642 a = self.alpha 

643 b = self.beta 

644 return a * b / ((a + b) ** 2 * (a + b + 1)) * self.length 

645 

646 @property 

647 def var(self) -> float: 

648 """Variance of the distribution.""" 

649 return self.cov 

650 

651 @property 

652 def std(self) -> float: 

653 """Standard deviation of the distribution.""" 

654 return np.sqrt(self.var) 

655 

656 def pdf(self, x: np.ndarray) -> np.ndarray: 

657 """Evaluate PDF. 

658 

659 Parameters 

660 ---------- 

661 x : array_like 

662 Evaluation points. 

663 

664 Returns 

665 ------- 

666 : 

667 Values of PDF evaluated in `x`. 

668 """ 

669 y = pt.misc.shift_coord(x, self.domain.flatten(), [0, 1]) 

670 ret = scipy.stats.beta.pdf(y, a=self.alpha, b=self.beta) 

671 return ret / self.length 

672 

673 def log_pdf(self, x: np.ndarray) -> np.ndarray: 

674 """Evaluate log-PDF. 

675 

676 Parameters 

677 ---------- 

678 x : array_like 

679 Evaluation points. 

680 

681 Returns 

682 ------- 

683 : 

684 Values of log-PDF evaluated in `x`. 

685 """ 

686 y = pt.misc.shift_coord(x, self.domain.flatten(), [0, 1]) 

687 ret = scipy.stats.beta.logpdf(y, a=self.alpha, b=self.beta) 

688 return ret - np.log(self.length) 

689 

690 def grad_x_log_pdf(self, x: np.ndarray) -> np.ndarray: 

691 """Evaluate gradient of log-PDF. 

692 

693 .. note:: 

694 Not yet implemented. 

695 

696 Parameters 

697 ---------- 

698 x : array_like 

699 Evaluation points. 

700 

701 Returns 

702 ------- 

703 : 

704 Values of gradient (vector valued) of log-PDF evaluated in `x`. 

705 """ 

706 raise NotImplementedError 

707 

708 def hess_x_log_pdf(self, x: np.ndarray) -> np.ndarray: 

709 """Evaluate Hessian of log-PDF. 

710 

711 .. note:: 

712 Not yet implemented. 

713 

714 Parameters 

715 ---------- 

716 x : array_like 

717 Evaluation points. 

718 

719 Returns 

720 ------- 

721 : 

722 Values of Hessian (matrix valued) of log-PDF evaluated in `x`. 

723 """ 

724 raise NotImplementedError 

725 

726 def sample(self, shape: int | list | tuple | np.ndarray) -> np.ndarray: 

727 """Draw samples from distribution. 

728 

729 Parameters 

730 ---------- 

731 shape : array_like 

732 Shape of the samples. 

733 

734 Returns 

735 ------- 

736 : 

737 Random samples of specified shape. 

738 """ 

739 SAMPLE = np.random.beta(self.alpha, self.beta, shape) 

740 return pt.misc.shift_coord(SAMPLE, [0, 1], self.domain.flatten()) 

741 

742 

743class WLSUnivariateSampler(Sampler): 

744 """Sampler for univariate optimally distributed samples on given domain. 

745 

746 Given a stochastic variable :math:`y\\in\\Gamma\\subset\\mathbb{R}` 

747 with :math:`y\\sim\\pi` and a finite subset 

748 :math:`\\{P_j\\}_{j=0}^{d-1}` of an orthonormal polynomial basis of 

749 :math:`L^2(\\Gamma,\\pi)`, the optimal weighted least-squares sampling 

750 distribution for a function 

751 :math:`u\\in\\operatorname{span}\\{P_j\\ \\vert\\ j=0,\\dots,d-1 \\}` reads 

752 

753 .. math:: 

754 \\mathrm{d}\\mu = w^{-1} \\mathrm{d}\\pi 

755 \\qquad\\mbox{with weight}\\qquad 

756 w^{-1}(y) = \\frac{1}{d}\\sum_{j=0}^{d-1}\\vert P_j(y)\\vert^2. 

757 

758 Parameters 

759 ---------- 

760 domain : array_like 

761 Interval of support of distribution. 

762 

763 Notes 

764 ----- 

765 To generate samples from the weighted least-squares distribution rejection 

766 sampling is used. For certain basis functions it is possible to choose a 

767 well-suited trial sampler for the rejection sampling, which can be enabled 

768 via setting ``tsa=True``. 

769 

770 See Also 

771 -------- 

772 pythia.sampler.WLSTensorSampler 

773 

774 References 

775 ---------- 

776 The optimal weighted least-squares sampling is based on the results in 

777 Cohen & Migliorati [1]_. 

778 

779 .. [1] Cohen, A. and Migliorati, G., 

780 “Optimal weighted least-squares methods”, 

781 SMAI Journal of Computational Mathematics 3, 181-203 (2017). 

782 """ 

783 

784 def __init__( 

785 self, param: pt.parameter.Parameter, deg: int, tsa: bool = True 

786 ) -> None: 

787 """Initiate WLSUnivariateSampler object.""" 

788 self.parameter = param 

789 self.deg = deg 

790 self._base_sampler = assign_sampler(param) 

791 self.domain = self._base_sampler.domain 

792 self._basis = pt.basis.univariate_basis([self.parameter], [self.deg])[0] 

793 self._tsa = tsa 

794 self._trial_sampler, self._bulk = self._compute_trial_sampler() 

795 assert self.domain.shape == (1, 2) 

796 # self._length = self.domain[0, 1] - self.domain[0, 0] 

797 

798 @property 

799 def dimension(self) -> int: 

800 """Parameter dimension.""" 

801 return self.domain.shape[0] 

802 

803 @property 

804 def mass(self) -> float: 

805 """Mass of the PDF.""" 

806 return 1 

807 

808 @property 

809 def maximum(self) -> float: 

810 """Maximum value of the PDF.""" 

811 if self.parameter.distribution == "uniform": 

812 # alternatively return pdf in domain[0, 1] 

813 return self.pdf(self.domain[0, 0]) 

814 elif self.parameter.distribution == "normal": 

815 if self.deg % 2 == 0: 

816 return self.pdf(self.parameter.mean) 

817 else: 

818 # Note: this can be improved 

819 return get_maximum(self.pdf, self.domain) 

820 elif self.parameter.distribution == "gamma": 

821 # Note: this is probably not going to work for domain = [a, inf] 

822 return get_maximum(self.pdf, self.domain) 

823 elif self.parameter.distribution == "beta": 

824 # Note: this can be improved 

825 return get_maximum(self.pdf, self.domain) 

826 return get_maximum(self.pdf, self.domain) 

827 

828 @property 

829 def mean(self) -> float: 

830 """Mean value of the distribution.""" 

831 mean, err = quad( 

832 lambda x: x * self.pdf(x), self.domain[0, 0], self.domain[0, 1] 

833 ) 

834 if err > 1e-8: 

835 warnings.warn(f"quadrature error large: ({err})") 

836 return mean 

837 

838 @property 

839 def cov(self) -> float: 

840 """(Co)Variance of the distribution.""" 

841 moment, err = quad( 

842 lambda x: x**2 * self.pdf(x), self.domain[0, 0], self.domain[0, 1] 

843 ) 

844 if err > 1e-8: 

845 warnings.warn(f"quadrature error large: ({err})") 

846 return moment - self.mean**2 

847 

848 @property 

849 def var(self) -> float: 

850 """Variance of the distribution.""" 

851 return self.cov 

852 

853 @property 

854 def std(self) -> float: 

855 """Standard deviation of the distribution.""" 

856 return np.sqrt(self.var) 

857 

858 def weight(self, x: np.ndarray | float | int) -> np.ndarray: 

859 """Weights for the pdf. 

860 

861 Parameters 

862 ---------- 

863 x : np.ndarray 

864 Points the weight function is evaluated in. 

865 

866 Returns 

867 ------- 

868 w : array_like 

869 Weights of evaluation points `x`. 

870 """ 

871 if isinstance(x, (float, int)): 

872 x = np.array(x).reshape(1) 

873 if x.ndim == 2 and x.shape[1] == 1: 

874 x = np.array(x).reshape(-1) 

875 assert x.ndim < 2 

876 b_eval = np.sum([np.abs(p(x)) ** 2 for p in self._basis], axis=0) 

877 return (self.deg + 1) / b_eval 

878 

879 def pdf(self, x: np.ndarray | float | int) -> np.ndarray: 

880 """Evaluate uniform PDF. 

881 

882 Parameters 

883 ---------- 

884 x : array_like 

885 Evaluation points. 

886 

887 Returns 

888 ------- 

889 : 

890 Values of PDF evaluated in `x`. 

891 """ 

892 if isinstance(x, (float, int)): 

893 x = np.array(x).reshape(1) 

894 if x.ndim == 2 and x.shape[1] == 1: 

895 x = np.array(x).reshape(-1) 

896 return self._base_sampler.pdf(x) / self.weight(x) 

897 

898 def log_pdf(self, x: np.ndarray | float | int) -> np.ndarray: 

899 """Evaluate uniform log-PDF. 

900 

901 Parameters 

902 ---------- 

903 x : array_like 

904 Evaluation points. 

905 

906 Returns 

907 ------- 

908 : 

909 Values of log-PDF evaluated in `x`. 

910 """ 

911 if isinstance(x, (float, int)): 

912 x = np.array(x).reshape(1) 

913 if x.ndim == 2 and x.shape[1] == 1: 

914 x = np.array(x).reshape(-1) 

915 return self._base_sampler.log_pdf(x) - np.log(self.weight(x)) 

916 

917 def grad_x_log_pdf(self, x: np.ndarray) -> np.ndarray: 

918 """Evaluate gradient of uniform log-PDF. 

919 

920 Parameters 

921 ---------- 

922 x : array_like 

923 Evaluation points. 

924 

925 Returns 

926 ------- 

927 : 

928 Values of gradient (vector valued) of log-PDF evaluated in `x`. 

929 """ 

930 raise NotImplementedError 

931 

932 def hess_x_log_pdf(self, x: np.ndarray) -> np.ndarray: 

933 """Evaluate Hessian of uniform log-PDF. 

934 

935 Parameters 

936 ---------- 

937 x : array_like 

938 Evaluation points. 

939 

940 Returns 

941 ------- 

942 : 

943 Values of Hessian (matrix valued) of log-PDF evaluated in `x`. 

944 """ 

945 raise NotImplementedError 

946 

947 def _compute_trial_sampler(self) -> tuple[Sampler, float]: 

948 """Trial sampler adaptation. 

949 

950 .. note:: 

951 TSA currently only available for uniform parameter distribution. 

952 

953 Parameters 

954 ---------- 

955 tsa : bool 

956 Adapt trial sampler or simply use uniform product sampler. 

957 

958 Returns 

959 ------- 

960 trial_sampler : pt.sampler.Sampler 

961 Trial sampler. 

962 bulk : array_like 

963 Domain estimate of the mass of the distribution. 

964 """ 

965 if self._tsa is True and self.parameter.distribution == "uniform": 

966 trial_sampler = BetaSampler(self.parameter.domain, 0.5, 0.5) 

967 else: 

968 trial_sampler = UniformSampler(self.parameter.domain) 

969 bulk = get_maximum(lambda x: self.pdf(x) / trial_sampler.pdf(x), self.domain) 

970 return trial_sampler, bulk 

971 

972 def sample(self, shape: int | list | tuple | np.ndarray) -> np.ndarray: 

973 """Draw samples from weighted least-squares parameter distribution. 

974 

975 Parameters 

976 ---------- 

977 shape : array_like 

978 Shape of the samples. 

979 

980 Returns 

981 ------- 

982 : 

983 Random samples of specified shape. 

984 """ 

985 if isinstance(shape, int): 

986 shape = (shape,) 

987 samples = rejection_sampling( 

988 self.pdf, self._trial_sampler, self._bulk, self.dimension, shape 

989 ) 

990 return samples 

991 

992 

993class ProductSampler(Sampler): 

994 """Tensor sampler for independent parameters. 

995 

996 Sampler for cartesian product samples of a list of (independent) univariate 

997 samplers. 

998 

999 Parameters 

1000 ---------- 

1001 sampler_list : list of `pythia.sampler.Sampler` 

1002 list of (univariate) Sampler objects. 

1003 """ 

1004 

1005 def __init__(self, sampler_list: list[Sampler]) -> None: 

1006 """Initiate ProductSampler object.""" 

1007 self.samplers = list(sampler_list) 

1008 self.domain = np.squeeze(np.array([s.domain for s in self.samplers])) 

1009 # Add dimension if len(sampler_list) == 1. 

1010 if self.domain.ndim < 2: 

1011 self.domain.shape = 1, 2 

1012 # sampler_list should contain 1D samplers only 

1013 assert self.domain.shape == (len(sampler_list), 2) 

1014 

1015 @property 

1016 def dimension(self) -> int: 

1017 """Dimension of the parameters.""" 

1018 return self.domain.shape[0] 

1019 

1020 @property 

1021 def mass(self) -> float: 

1022 """Mass of the PDF.""" 

1023 return np.prod(np.array([s.mass for s in self.samplers])) 

1024 

1025 @property 

1026 def maximum(self) -> float: 

1027 """Maximum value of the PDF.""" 

1028 return np.prod(np.array([s.maximum for s in self.samplers])) 

1029 

1030 @property 

1031 def mean(self) -> np.ndarray: 

1032 """Mean of the PDF.""" 

1033 return np.array([s.mean for s in self.samplers]) 

1034 

1035 @property 

1036 def cov(self) -> np.ndarray: 

1037 """Covariance of the PDF.""" 

1038 return np.diag([s.cov for s in self.samplers]) 

1039 

1040 def weight(self, x: np.ndarray) -> np.ndarray: 

1041 """Weights of the product PDF. 

1042 

1043 Parameters 

1044 ---------- 

1045 x : np.ndarray 

1046 Evaluation points. 

1047 

1048 Returns 

1049 ------- 

1050 : 

1051 Array of uniform weights for samples. 

1052 """ 

1053 if x.ndim == 1: 

1054 x.shape = 1, -1 

1055 assert x.ndim == 2 

1056 return np.ones(x.shape[0]) 

1057 

1058 def pdf(self, x: np.ndarray) -> np.ndarray: 

1059 """Evaluate PDF. 

1060 

1061 The PDF is given by the product of the univariate PDFs. 

1062 

1063 Parameters 

1064 ---------- 

1065 x : array_like 

1066 Evaluation points. 

1067 

1068 Returns 

1069 ------- 

1070 : 

1071 Values of PDF evaluated in `x`. 

1072 """ 

1073 assert x.shape[-1] == self.dimension 

1074 densities = [s.pdf for s in self.samplers] 

1075 val = np.array([densities[jj](x[..., jj]) for jj in range(self.dimension)]) 

1076 return np.prod(val, axis=0) 

1077 

1078 def log_pdf(self, x: np.ndarray) -> np.ndarray: 

1079 """Evaluate log-PDF. 

1080 

1081 The log-PDF is given by the sum of the univariate log-PDFs. 

1082 

1083 Parameters 

1084 ---------- 

1085 x : array_like 

1086 Evaluation points. 

1087 

1088 Returns 

1089 ------- 

1090 : 

1091 Values of log-PDF evaluated in `x`. 

1092 """ 

1093 assert x.shape[-1] == self.dimension 

1094 densities = [s.log_pdf for s in self.samplers] 

1095 val = np.array([densities[jj](x[..., jj]) for jj in range(self.dimension)]) 

1096 return np.sum(val, axis=0) 

1097 

1098 def grad_x_log_pdf(self, x: np.ndarray) -> np.ndarray: 

1099 """Evaluate gradient of log-PDF. 

1100 

1101 Parameters 

1102 ---------- 

1103 x : array_like 

1104 Evaluation points. 

1105 

1106 Returns 

1107 ------- 

1108 : 

1109 Values of gradient (vector valued) of log-PDF evaluated in `x`. 

1110 """ 

1111 grad_densities = [s.grad_x_log_pdf for s in self.samplers] 

1112 return np.array( 

1113 [grad_densities[jj](x[..., jj]) for jj in range(self.dimension)] 

1114 ).T 

1115 

1116 def hess_x_log_pdf(self, x: np.ndarray) -> np.ndarray: 

1117 """Evaluate Hessian of log-PDF. 

1118 

1119 Parameters 

1120 ---------- 

1121 x : array_like 

1122 Evaluation points. 

1123 

1124 Returns 

1125 ------- 

1126 : 

1127 Values of Hessian (matrix valued) of log-PDF evaluated in `x`. 

1128 """ 

1129 # create an (x.shape[0], self._dim, self._dim) tensor where each 

1130 # (self._dim, self._dim) matrix is the identity 

1131 eye = np.tile(np.expand_dims(np.eye(self.dimension), 0), (x.shape[0], 1, 1)) 

1132 

1133 # create an (x.shape[0],self._dim,1) tensor where each (x.shape[0],1,1) 

1134 # subtensor is a diagonal entry of the hessian 

1135 hess_densities = [s.hess_x_log_pdf for s in self.samplers] 

1136 hess = np.expand_dims( 

1137 np.array( 

1138 [hess_densities[jj](x[..., jj]) for jj in range(self.dimension)] 

1139 ).T, 

1140 2, 

1141 ) 

1142 return eye * hess 

1143 

1144 def sample(self, shape: int | list | tuple | np.ndarray) -> np.ndarray: 

1145 """Draw samples from distribution. 

1146 

1147 Parameters 

1148 ---------- 

1149 shape : array_like 

1150 Shape of the samples. 

1151 

1152 Returns 

1153 ------- 

1154 : 

1155 Random samples of specified shape. 

1156 """ 

1157 if isinstance(shape, int): 

1158 shape = (shape,) 

1159 samples = [s.sample(shape) for s in self.samplers] 

1160 return np.stack(samples, -1) 

1161 

1162 

1163class ParameterSampler(Sampler): 

1164 """Product sampler of given parameters. 

1165 

1166 Parameters 

1167 ---------- 

1168 params : list of `pythia.parameter.Parameter` 

1169 list containing information of parameters. 

1170 """ 

1171 

1172 def __init__(self, params: list[pt.parameter.Parameter]) -> None: 

1173 """Initiate ParameterSampler object.""" 

1174 self.parameter = params 

1175 assert isinstance(self.parameter, list) 

1176 self.domain = np.array([param.domain for param in self.parameter]) 

1177 self._product_sampler = ProductSampler( 

1178 [assign_sampler(param) for param in self.parameter] 

1179 ) 

1180 

1181 @property 

1182 def dimension(self) -> int: 

1183 """Dimension of the parameters.""" 

1184 return self._product_sampler.dimension 

1185 

1186 @property 

1187 def mass(self) -> float: 

1188 """Mass of the PDF.""" 

1189 return self._product_sampler.mass 

1190 

1191 @property 

1192 def maximum(self) -> float: 

1193 """Maximum value of the PDF.""" 

1194 return self._product_sampler.maximum 

1195 

1196 @property 

1197 def mean(self) -> np.ndarray: 

1198 """Mean of the PDF.""" 

1199 return self._product_sampler.mean 

1200 

1201 @property 

1202 def cov(self) -> np.ndarray: 

1203 """Covariance of the PDF.""" 

1204 return self._product_sampler.cov 

1205 

1206 def weight(self, x: np.ndarray) -> np.ndarray: 

1207 """Weights of the parameter product PDF. 

1208 

1209 Parameters 

1210 ---------- 

1211 x : np.ndarray 

1212 Evaluation points. 

1213 

1214 Returns 

1215 ------- 

1216 : 

1217 Array of uniform weights for samples. 

1218 """ 

1219 return self._product_sampler.weight(x) 

1220 

1221 def pdf(self, x: np.ndarray) -> np.ndarray: 

1222 """Evaluate PDF. 

1223 

1224 Parameters 

1225 ---------- 

1226 x : array_like 

1227 Evaluation points. 

1228 

1229 Returns 

1230 ------- 

1231 : 

1232 Values of PDF evaluated in `x`. 

1233 """ 

1234 return self._product_sampler.pdf(x) 

1235 

1236 def log_pdf(self, x: np.ndarray) -> np.ndarray: 

1237 """Evaluate log-PDF. 

1238 

1239 The log-PDF is given by the sum of the univariate log-PDFs. 

1240 

1241 Parameters 

1242 ---------- 

1243 x : array_like 

1244 Evaluation points. 

1245 

1246 Returns 

1247 ------- 

1248 : 

1249 Values of log-PDF evaluated in `x`. 

1250 """ 

1251 return self._product_sampler.log_pdf(x) 

1252 

1253 def grad_x_log_pdf(self, x: np.ndarray) -> np.ndarray: 

1254 """Evaluate gradient of log-PDF. 

1255 

1256 Parameters 

1257 ---------- 

1258 x : array_like 

1259 Evaluation points. 

1260 

1261 Returns 

1262 ------- 

1263 : 

1264 Values of gradient (vector valued) of log-PDF evaluated in `x`. 

1265 """ 

1266 return self._product_sampler.grad_x_log_pdf(x) 

1267 

1268 def hess_x_log_pdf(self, x: np.ndarray) -> np.ndarray: 

1269 """Evaluate Hessian of log-PDF. 

1270 

1271 Parameters 

1272 ---------- 

1273 x : array_like 

1274 Evaluation points. 

1275 

1276 Returns 

1277 ------- 

1278 : 

1279 Values of Hessian (matrix valued) of log-PDF evaluated in `x`. 

1280 """ 

1281 return self._product_sampler.hess_x_log_pdf(x) 

1282 

1283 def sample(self, shape: int | list | tuple | np.ndarray) -> np.ndarray: 

1284 """Draw samples from distribution. 

1285 

1286 Parameters 

1287 ---------- 

1288 shape : array_like 

1289 Shape of the samples. 

1290 

1291 Returns 

1292 ------- 

1293 : 

1294 Random samples of specified shape. 

1295 """ 

1296 return self._product_sampler.sample(shape) 

1297 

1298 

1299class WLSSampler(Sampler): 

1300 """Weighted Least-Squares sampler. 

1301 

1302 Given a stochastic variable :math:`y\\in\\Gamma\\subset\\mathbb{R}^M` 

1303 with :math:`y\\sim\\pi`, a set of multiindices 

1304 :math:`\\Lambda\\subset\\mathbb{N}_0^M` and a finite subset 

1305 :math:`\\{P_\\alpha\\}_{\\alpha\\in\\Lambda}` of an orthonormal polynomial 

1306 basis of :math:`L^2(\\Gamma,\\pi)`, the optimal weighted least-squares 

1307 sampling distribution for a function 

1308 :math:`u\\in\\operatorname{span}\\{P_\\alpha\\ \\vert\\ \\alpha\\in\\Lambda\\}` 

1309 reads 

1310 

1311 .. math:: 

1312 \\mathrm{d}\\mu = w^{-1} \\mathrm{d}\\pi 

1313 \\qquad\\mbox{with weight}\\qquad 

1314 w^{-1}(y) = \\frac{1}{\\vert\\Lambda\\vert}\\sum_{\\alpha\\in\\Lambda}\\vert P_\\alpha(y)\\vert^2, 

1315 

1316 where :math:`\\vert\\Lambda\\vert` denotes the number of elements in 

1317 :math:`\\Lambda`. 

1318 

1319 Parameters 

1320 ---------- 

1321 params : list of `pythia.parameter.Parameter` 

1322 list of parameters. 

1323 basis : list 

1324 list of basis functions. 

1325 tsa : bool, default=False 

1326 Trial sampler adaptation. If True, a trial sampler is chosen on the 

1327 distributions of parameters, if false a uniform trial sampler is used. 

1328 

1329 Other Parameters 

1330 ---------------- 

1331 trial_sampler : pythia.sampler.Sampler, default=None 

1332 Trial sampler for rejection sampling. If `tsa` is true and either 

1333 `trial_sampler` or `bulk` are `None`, the trial sampler is chosen 

1334 automatically. 

1335 bulk : float, defaul=None 

1336 Scaling for trial sampler. If `tsa` is true and either 

1337 `trial_sampler` or `bulk` are `None`, the trial sampler is chosen 

1338 automatically. 

1339 

1340 Notes 

1341 ----- 

1342 To generate samples from the weighted least-squares distribution rejection 

1343 sampling is used. For certain basis functions it is possible to choose a 

1344 well-suited trial sampler for the rejection sampling, which can be enabled 

1345 via setting ``tsa=True``. 

1346 

1347 See Also 

1348 -------- 

1349 pythia.sampler.WLSUnivariateSampler, pythia.sampler.WLSTensorSampler 

1350 

1351 References 

1352 ---------- 

1353 The optimal weighted least-squares sampling is based on the results of 

1354 Cohen & Migliorati [1]_. 

1355 """ 

1356 

1357 def __init__( 

1358 self, 

1359 params: list[pt.parameter.Parameter], 

1360 basis: list[Callable], 

1361 tsa: bool = True, 

1362 trial_sampler: Sampler | None = None, 

1363 bulk: float | None = None, 

1364 ) -> None: 

1365 """Initiate WLSSampler object.""" 

1366 self.parameter = params 

1367 self._param_sampler = ParameterSampler(self.parameter) 

1368 self.domain = self._param_sampler.domain 

1369 self.basis = basis 

1370 self._tsa = tsa 

1371 if trial_sampler is None or bulk is None: 

1372 self._trial_sampler, self._bulk = self._compute_trial_sampler() 

1373 else: 

1374 self._trial_sampler = trial_sampler 

1375 self._bulk = bulk 

1376 

1377 @property 

1378 def dimension(self) -> int: 

1379 """Dimension of the parameters.""" 

1380 return self._param_sampler.dimension 

1381 

1382 @property 

1383 def mass(self): 

1384 """Mass of the PDF.""" 

1385 return 1 

1386 

1387 @property 

1388 def maximum(self) -> float: 

1389 """Maximum value of the PDF.""" 

1390 return get_maximum(self.pdf, self.domain) 

1391 

1392 @property 

1393 def mean(self) -> np.ndarray: 

1394 """Mean of the PDF.""" 

1395 raise NotImplementedError 

1396 

1397 @property 

1398 def cov(self) -> np.ndarray: 

1399 """Covariance of the PDF.""" 

1400 raise NotImplementedError 

1401 

1402 def weight(self, x: np.ndarray) -> np.ndarray: 

1403 """Weights for the PDF. 

1404 

1405 Parameters 

1406 ---------- 

1407 x : array_like 

1408 Points the weight function is evaluated in. 

1409 

1410 Returns 

1411 ------- 

1412 : 

1413 weights of evaluation points `x`. 

1414 """ 

1415 km = np.sum([abs(p(x)) ** 2 for p in self.basis], axis=0) 

1416 return len(self.basis) / km 

1417 

1418 def pdf(self, x: np.ndarray) -> np.ndarray: 

1419 """Evaluate PDF. 

1420 

1421 Parameters 

1422 ---------- 

1423 x : array_like 

1424 Evaluation points. 

1425 

1426 Returns 

1427 ------- 

1428 : 

1429 Values of PDF evaluated in `x`. 

1430 """ 

1431 assert x.shape[-1] == self.dimension 

1432 return self._param_sampler.pdf(x) / self.weight(x) 

1433 

1434 def log_pdf(self, x: np.ndarray) -> np.ndarray: 

1435 """Evaluate log-PDF. 

1436 

1437 The log-PDF is given by the sum of the univariate log-PDFs. 

1438 

1439 Parameters 

1440 ---------- 

1441 x : array_like 

1442 Evaluation points. 

1443 

1444 Returns 

1445 ------- 

1446 : 

1447 Values of log-PDF evaluated in `x`. 

1448 """ 

1449 assert x.shape[-1] == self.dimension 

1450 return self._param_sampler.log_pdf(x) - np.log(self.weight(x)) 

1451 

1452 def grad_x_log_pdf(self, x: np.ndarray) -> np.ndarray: 

1453 """Evaluate gradient of log-PDF. 

1454 

1455 Parameters 

1456 ---------- 

1457 x : array_like 

1458 Evaluation points. 

1459 

1460 Returns 

1461 ------- 

1462 : 

1463 Values of gradient (vector valued) of log-PDF evaluated in `x`. 

1464 """ 

1465 raise NotImplementedError 

1466 

1467 def hess_x_log_pdf(self, x: np.ndarray) -> np.ndarray: 

1468 """Evaluate Hessian of log-PDF. 

1469 

1470 Parameters 

1471 ---------- 

1472 x : array_like 

1473 Evaluation points. 

1474 

1475 Returns 

1476 ------- 

1477 : 

1478 Values of Hessian (matrix valued) of log-PDF evaluated in `x`. 

1479 """ 

1480 raise NotImplementedError 

1481 

1482 def _compute_trial_sampler(self) -> tuple[Sampler, float]: 

1483 """Trial sampler adaptation. 

1484 

1485 Parameters 

1486 ---------- 

1487 tsa : bool 

1488 Adapt trial sampler or simply use uniform product sampler. 

1489 

1490 Returns 

1491 ------- 

1492 trial_sampler : pt.sampler.Sampler 

1493 Trial sampler. 

1494 bulk : array_like 

1495 Domain estimate of the mass of the distribution. 

1496 """ 

1497 sampler_list = [] 

1498 for param in self.parameter: 

1499 if self._tsa is True and param.distribution == "uniform": 

1500 sampler_list.append(BetaSampler(param.domain, 0.5, 0.5)) 

1501 else: 

1502 sampler_list.append(UniformSampler(param.domain)) 

1503 trial_sampler = ProductSampler(sampler_list) 

1504 bulk = get_maximum(lambda x: self.pdf(x) / trial_sampler.pdf(x), self.domain) 

1505 return trial_sampler, bulk 

1506 

1507 def sample(self, shape: int | list | tuple | np.ndarray) -> np.ndarray: 

1508 """Draw samples from distribution. 

1509 

1510 Parameters 

1511 ---------- 

1512 shape : array_like 

1513 Shape of the samples. 

1514 

1515 Returns 

1516 ------- 

1517 : 

1518 Random samples of specified shape. 

1519 """ 

1520 if isinstance(shape, int): 

1521 shape = (shape,) 

1522 samples = rejection_sampling( 

1523 self.pdf, self._trial_sampler, self._bulk, self.dimension, shape 

1524 ) 

1525 return samples 

1526 

1527 

1528class WLSTensorSampler(Sampler): 

1529 """Weighted least-squares sampler for tensor multivariate basis. 

1530 

1531 Given a stochastic variable :math:`y\\in\\Gamma\\subset\\mathbb{R}^M` 

1532 with :math:`y\\sim\\pi=\\prod_{m=1}^M\\pi_m` for one dimensional densities 

1533 :math:`\\pi_m`, a tensor set of multiindices 

1534 :math:`\\Lambda=[d_1]\\times\\dots\\times[d_M]\\subset\\mathbb{N}_0^M`, 

1535 where :math:`[d_m]=\\{0,\\dots,d_m-1\\}`, and a finite subset 

1536 :math:`\\{P_\\alpha\\}_{\\alpha\\in\\Lambda}` of an orthonormal product 

1537 polynomial basis of :math:`L^2(\\Gamma,\\pi)`, i.e., 

1538 :math:`P_\\alpha(y) = \\prod_{m=1}^M P_{\\alpha_m}(y_m)`, the optimal 

1539 weighted least-squares sampling distribution for a function 

1540 :math:`u\\in\\operatorname{span}\\{P_\\alpha\\ \\vert\\ \\alpha\\in\\Lambda\\}` 

1541 reads 

1542 

1543 .. math:: 

1544 \\mathrm{d}\\mu = w^{-1} \\mathrm{d}\\pi 

1545 \\qquad\\mbox{with weight}\\qquad 

1546 w^{-1}(y) = \\prod_{m=1}^M\\frac{1}{d_m}\\sum_{\\alpha_m\\in[d_m]}\\vert P_{\\alpha_m}(y_m)\\vert^2. 

1547 

1548 Parameters 

1549 ---------- 

1550 params : list of `pythia.parameter.Parameter` 

1551 Parameter list. 

1552 deg : list of int 

1553 Polynomial degree of each component (same for all). 

1554 tsa : bool, default=True 

1555 Trial sampler adaptation. If True, a trial sampler is chosen on the 

1556 distributions of parameters, if false a uniform trial sampler is 

1557 used. 

1558 

1559 Notes 

1560 ----- 

1561 To generate samples from the weighted least-squares distribution rejection 

1562 sampling is used. For certain basis functions it is possible to choose a 

1563 well-suited trial sampler for the rejection sampling, which can be enabled 

1564 via setting ``tsa=True``. 

1565 

1566 See Also 

1567 -------- 

1568 pythia.sampler.WLSUnivariateSampler 

1569 

1570 References 

1571 ---------- 

1572 The optimal weighted least-squares sampling is based on the results in 

1573 Cohen & Migliorati [1]_. 

1574 """ 

1575 

1576 def __init__( 

1577 self, params: list[pt.parameter.Parameter], deg: list[int], tsa: bool = True 

1578 ) -> None: 

1579 """Initiate WLSTensorSampler object.""" 

1580 self.parameter = params 

1581 self.deg = deg 

1582 self.domain = np.array([param.domain for param in self.parameter]) 

1583 self._tsa = tsa 

1584 # self._uBasis = pt.basis.univariate_basis(self.parameter, self.deg) 

1585 # self._univariate_samplers = [assign_sampler(param) 

1586 # for param in self.parameter] 

1587 self._univariate_samplers = [ 

1588 WLSUnivariateSampler(p, d, self._tsa) 

1589 for (p, d) in zip(self.parameter, self.deg) 

1590 ] 

1591 self._product_sampler = ProductSampler(self._univariate_samplers) 

1592 

1593 @property 

1594 def dimension(self) -> int: 

1595 """Dimension of the parameters.""" 

1596 return self._product_sampler.dimension 

1597 

1598 @property 

1599 def mass(self) -> float: 

1600 """Mass of the PDF.""" 

1601 return self._product_sampler.mass 

1602 

1603 @property 

1604 def maximum(self) -> float: 

1605 """Maximum value of the PDF.""" 

1606 return self._product_sampler.maximum 

1607 

1608 @property 

1609 def mean(self) -> np.ndarray: 

1610 """Mean of the PDF.""" 

1611 return self._product_sampler.mean 

1612 

1613 @property 

1614 def cov(self) -> np.ndarray: 

1615 """Covariance of the PDF.""" 

1616 return self._product_sampler.cov 

1617 

1618 def weight(self, x: np.ndarray) -> np.ndarray: 

1619 """Weights for the PDF. 

1620 

1621 Parameters 

1622 ---------- 

1623 x : array_like 

1624 Points the weight function is evaluated in. 

1625 

1626 Returns 

1627 ------- 

1628 : 

1629 Weights of evaluation points `x`. 

1630 """ 

1631 assert x.shape[-1] == self.dimension 

1632 weights = [s.weight for s in self._univariate_samplers] 

1633 val = np.array([weights[jj](x[..., jj]) for jj in range(self.dimension)]) 

1634 return np.prod(val, axis=0) 

1635 

1636 def pdf(self, x: np.ndarray) -> np.ndarray: 

1637 """Evaluate PDF. 

1638 

1639 Parameters 

1640 ---------- 

1641 x : array_like 

1642 Evaluation points. 

1643 

1644 Returns 

1645 ------- 

1646 : 

1647 Values of PDF evaluated in `x`. 

1648 """ 

1649 return self._product_sampler.pdf(x) 

1650 

1651 def log_pdf(self, x: np.ndarray) -> np.ndarray: 

1652 """Evaluate log-PDF. 

1653 

1654 The log-PDF is given by the sum of the univariate log-PDFs. 

1655 

1656 Parameters 

1657 ---------- 

1658 x : array_like 

1659 Evaluation points. 

1660 

1661 Returns 

1662 ------- 

1663 : 

1664 Values of log-PDF evaluated in `x`. 

1665 """ 

1666 return self._product_sampler.log_pdf(x) 

1667 

1668 def grad_x_log_pdf(self, x: np.ndarray) -> np.ndarray: 

1669 """Evaluate gradient of log-PDF. 

1670 

1671 Parameters 

1672 ---------- 

1673 x : array_like 

1674 Evaluation points. 

1675 

1676 Returns 

1677 ------- 

1678 : 

1679 Values of gradient (vector valued) of log-PDF evaluated in `x`. 

1680 """ 

1681 return self._product_sampler.grad_x_log_pdf(x) 

1682 

1683 def hess_x_log_pdf(self, x: np.ndarray) -> np.ndarray: 

1684 """Evaluate Hessian of log-PDF. 

1685 

1686 Parameters 

1687 ---------- 

1688 x : array_like 

1689 Evaluation points. 

1690 

1691 Returns 

1692 ------- 

1693 : 

1694 Values of Hessian (matrix valued) of log-PDF evaluated in `x`. 

1695 """ 

1696 return self._product_sampler.hess_x_log_pdf(x) 

1697 

1698 def sample(self, shape: int | list | tuple | np.ndarray) -> np.ndarray: 

1699 """Draw samples from distribution. 

1700 

1701 Parameters 

1702 ---------- 

1703 shape : array_like 

1704 Shape of the samples. 

1705 

1706 Returns 

1707 ------- 

1708 : 

1709 Random samples of specified shape. 

1710 """ 

1711 return self._product_sampler.sample(shape) 

1712 

1713 

1714def rejection_sampling( 

1715 pdf: Callable, 

1716 trial_sampler: Sampler, 

1717 scale: float, 

1718 dimension: int, 

1719 shape: int | list | tuple | np.ndarray, 

1720) -> np.ndarray: 

1721 """Draw samples from pdf by rejection sampling. 

1722 

1723 Parameters 

1724 ---------- 

1725 pdf : Callable 

1726 Probability density the samples are generated from. 

1727 trial_sampler : Sampler 

1728 Trial sampler proposal samples are drawn from. 

1729 scale : float 

1730 Threshold parameter with ``pdf <= scale * trialSampler.pdf`` 

1731 dimension : int 

1732 Dimension of the (input of the) pdf. 

1733 shape : array_like 

1734 Shape of the samples. 

1735 

1736 Returns 

1737 ------- 

1738 : 

1739 Random samples of specified shape. 

1740 """ 

1741 if isinstance(shape, int): 

1742 shape = (shape,) 

1743 size = np.prod(np.array(shape)) 

1744 trial_samples = trial_sampler.sample(size) 

1745 samples = np.empty_like(trial_samples) 

1746 pointer = 0 

1747 while pointer < size: 

1748 trial_samples = trial_sampler.sample(size) 

1749 checker = np.random.rand(trial_samples.shape[0]) 

1750 is_valid = checker * max(scale, 1) * trial_sampler.pdf(trial_samples) <= pdf( 

1751 trial_samples 

1752 ) 

1753 valid_samples = trial_samples[is_valid] 

1754 start = pointer 

1755 end = min(pointer + valid_samples.shape[0], size) 

1756 samples[start:end] = valid_samples[ 

1757 0 : min(valid_samples.shape[0], size - pointer) 

1758 ] 

1759 pointer = end 

1760 if dimension > 1: 

1761 return np.moveaxis(samples.T.reshape(dimension, *shape), 0, -1) 

1762 return samples.reshape(*shape) 

1763 

1764 

1765def constraint_sampling( 

1766 sampler: Sampler, 

1767 constraints: list[Callable], 

1768 shape: int | list | tuple | np.ndarray, 

1769) -> np.ndarray: 

1770 """Draw samples according to algebraic constraints. 

1771 

1772 Draw samples from target distribution and discard samples that do not 

1773 satisfy the constraints. 

1774 

1775 Parameters 

1776 ---------- 

1777 sampler : Sampler 

1778 Sampler to sample from. 

1779 constraints : list of callable 

1780 list of functions that return True if sample point satisfies the 

1781 constraint. 

1782 

1783 Returns 

1784 ------- 

1785 : 

1786 Samples drawn from sampler satisfying the constraints. 

1787 

1788 Notes 

1789 ----- 

1790 The constaints may lead to a non-normalized density function. 

1791 """ 

1792 if isinstance(shape, int): 

1793 shape = (shape,) 

1794 samples = np.empty(sampler.sample(np.prod(shape)).shape) 

1795 for j_sample in range(samples.shape[0]): 

1796 proposal = sampler.sample((1,)) 

1797 met_constraints = np.all([c(proposal) for c in constraints]) 

1798 while not met_constraints: 

1799 proposal = sampler.sample((1,)) 

1800 met_constraints = np.all([c(proposal) for c in constraints]) 

1801 samples[j_sample] = proposal 

1802 return np.moveaxis(samples.T.reshape(sampler.dimension, *shape), 0, -1) 

1803 

1804 

1805def assign_sampler(param: pt.parameter.Parameter) -> Sampler: 

1806 """Assign a univariate sampler to the given parameter. 

1807 

1808 Parameters 

1809 ---------- 

1810 param : pythia.parameter.Parameter 

1811 

1812 Returns 

1813 ------- 

1814 : 

1815 Univariate sampler. 

1816 """ 

1817 if param.distribution == "uniform": 

1818 return UniformSampler(param.domain) 

1819 elif param.distribution == "normal": 

1820 return NormalSampler(param.mean, param.var) 

1821 elif param.distribution == "gamma": 

1822 return GammaSampler(param.domain, param.alpha, param.beta) 

1823 elif param.distribution == "beta": 

1824 return BetaSampler(param.domain, param.alpha, param.beta) 

1825 raise ValueError(f"unknown distribution '{param.distribution}'") 

1826 

1827 

1828def get_maximum(f: Callable, domain: list | tuple | np.ndarray, n: int = 1000) -> float: 

1829 """Compute essential maximum of function by point evaluations. 

1830 

1831 Parameters 

1832 ---------- 

1833 f : callable 

1834 Function to evaluate. Needs to map from n-dim space to 1-dim space. 

1835 domain : array_like 

1836 Domain to evaluate function on. 

1837 n : int, default=1000 

1838 Number of function evaluations. Evaluations are done on a uniform grid 

1839 in domain. Actual number of points may thus be a little greater. 

1840 

1841 Returns 

1842 ------- 

1843 : 

1844 Approximation of maximum of function `f`. 

1845 """ 

1846 if not isinstance(domain, np.ndarray): 

1847 domain = np.array(domain) 

1848 assert domain.shape[-1] == 2 and domain.ndim == 2 

1849 n_points = int(np.ceil(np.power(n, 1 / domain.shape[0]))) 

1850 if domain.shape[0] > 1: 

1851 eps = np.finfo(float).eps # circumvent inf on bdry 

1852 x = pt.misc.cart_prod( 

1853 [np.linspace(dom[0] + eps, dom[1] - eps, n_points) for dom in domain] 

1854 ) 

1855 else: 

1856 x = np.linspace(domain[0, 0], domain[0, 1], n_points).reshape(-1, 1) 

1857 return np.max(np.abs(f(x))) 

1858 

1859 

1860if __name__ == "__main__": 

1861 pass