1from collections.abc import MutableSequence 2import warnings 3import io 4import copy 5 6import numpy as np 7import pandas as pd 8 9from . import endf 10import openmc.checkvalue as cv 11from .resonance import Resonances 12 13 14def _add_file2_contributions(file32params, file2params): 15 """Function for aiding in adding resonance parameters from File 2 that are 16 not always present in File 32. Uses already imported resonance data. 17 18 Paramaters 19 ---------- 20 file32params : pandas.Dataframe 21 Incomplete set of resonance parameters contained in File 32. 22 file2params : pandas.Dataframe 23 Resonance parameters from File 2. Ordered by energy. 24 25 Returns 26 ------- 27 parameters : pandas.Dataframe 28 Complete set of parameters ordered by L-values and then energy 29 30 """ 31 # Use l-values and competitiveWidth from File 2 data 32 # Re-sort File 2 by energy to match File 32 33 file2params = file2params.sort_values(by=['energy']) 34 file2params.reset_index(drop=True, inplace=True) 35 # Sort File 32 parameters by energy as well (maintaining index) 36 file32params.sort_values(by=['energy'], inplace=True) 37 # Add in values (.values converts to array first to ignore index) 38 file32params['L'] = file2params['L'].values 39 if 'competitiveWidth' in file2params.columns: 40 file32params['competitiveWidth'] = file2params['competitiveWidth'].values 41 # Resort to File 32 order (by L then by E) for use with covariance 42 file32params.sort_index(inplace=True) 43 return file32params 44 45 46class ResonanceCovariances(Resonances): 47 """Resolved resonance covariance data 48 49 Parameters 50 ---------- 51 ranges : list of openmc.data.ResonanceCovarianceRange 52 Distinct energy ranges for resonance data 53 54 Attributes 55 ---------- 56 ranges : list of openmc.data.ResonanceCovarianceRange 57 Distinct energy ranges for resonance data 58 59 """ 60 61 @property 62 def ranges(self): 63 return self._ranges 64 65 @ranges.setter 66 def ranges(self, ranges): 67 cv.check_type('resonance ranges', ranges, MutableSequence) 68 self._ranges = cv.CheckedList(ResonanceCovarianceRange, 69 'resonance range', ranges) 70 71 @classmethod 72 def from_endf(cls, ev, resonances): 73 """Generate resonance covariance data from an ENDF evaluation. 74 75 Parameters 76 ---------- 77 ev : openmc.data.endf.Evaluation 78 ENDF evaluation 79 resonances : openmc.data.Resonance object 80 openmc.data.Resonanance object generated from the same evaluation 81 used to import values not contained in File 32 82 83 Returns 84 ------- 85 openmc.data.ResonanceCovariances 86 Resonance covariance data 87 88 """ 89 file_obj = io.StringIO(ev.section[32, 151]) 90 91 # Determine whether discrete or continuous representation 92 items = endf.get_head_record(file_obj) 93 n_isotope = items[4] # Number of isotopes 94 95 ranges = [] 96 for iso in range(n_isotope): 97 items = endf.get_cont_record(file_obj) 98 abundance = items[1] 99 fission_widths = (items[3] == 1) # Flag for fission widths 100 n_ranges = items[4] # Number of resonance energy ranges 101 102 for j in range(n_ranges): 103 items = endf.get_cont_record(file_obj) 104 # Unresolved flags - 0: only scattering radius given 105 # 1: resolved parameters given 106 # 2: unresolved parameters given 107 unresolved_flag = items[2] 108 formalism = items[3] # resonance formalism 109 110 # Throw error for unsupported formalisms 111 if formalism in [0, 7]: 112 error = 'LRF='+str(formalism)+' covariance not supported '\ 113 'for this formalism' 114 raise NotImplementedError(error) 115 116 if unresolved_flag in (0, 1): 117 # Resolved resonance region 118 resonance = resonances.ranges[j] 119 erange = _FORMALISMS[formalism].from_endf(ev, file_obj, 120 items, resonance) 121 ranges.append(erange) 122 123 elif unresolved_flag == 2: 124 warn = 'Unresolved resonance not supported. Covariance '\ 125 'values for the unresolved region not imported.' 126 warnings.warn(warn) 127 128 return cls(ranges) 129 130 131class ResonanceCovarianceRange: 132 """Resonace covariance range. Base class for different formalisms. 133 134 Parameters 135 ---------- 136 energy_min : float 137 Minimum energy of the resolved resonance range in eV 138 energy_max : float 139 Maximum energy of the resolved resonance range in eV 140 141 Attributes 142 ---------- 143 energy_min : float 144 Minimum energy of the resolved resonance range in eV 145 energy_max : float 146 Maximum energy of the resolved resonance range in eV 147 parameters : pandas.DataFrame 148 Resonance parameters 149 covariance : numpy.array 150 The covariance matrix contained within the ENDF evaluation 151 lcomp : int 152 Flag indicating format of the covariance matrix within the ENDF file 153 file2res : openmc.data.ResonanceRange object 154 Corresponding resonance range with File 2 data. 155 mpar : int 156 Number of parameters in covariance matrix for each individual resonance 157 formalism : str 158 String descriptor of formalism 159 """ 160 def __init__(self, energy_min, energy_max): 161 self.energy_min = energy_min 162 self.energy_max = energy_max 163 164 def subset(self, parameter_str, bounds): 165 """Produce a subset of resonance parameters and the corresponding 166 covariance matrix to an IncidentNeutron object. 167 168 Parameters 169 ---------- 170 parameter_str : str 171 parameter to be discriminated 172 (i.e. 'energy', 'captureWidth', 'fissionWidthA'...) 173 bounds : np.array 174 [low numerical bound, high numerical bound] 175 176 Returns 177 ------- 178 res_cov_range : openmc.data.ResonanceCovarianceRange 179 ResonanceCovarianceRange object that contains a subset of the 180 covariance matrix (upper triangular) as well as a subset parameters 181 within self.file2params 182 183 """ 184 # Copy range and prevent change of original 185 res_cov_range = copy.deepcopy(self) 186 187 parameters = self.file2res.parameters 188 cov = res_cov_range.covariance 189 mpar = res_cov_range.mpar 190 # Create mask 191 mask1 = parameters[parameter_str] >= bounds[0] 192 mask2 = parameters[parameter_str] <= bounds[1] 193 mask = mask1 & mask2 194 res_cov_range.parameters = parameters[mask] 195 indices = res_cov_range.parameters.index.values 196 # Build subset of covariance 197 sub_cov_dim = len(indices)*mpar 198 cov_subset_vals = [] 199 for index1 in indices: 200 for i in range(mpar): 201 for index2 in indices: 202 for j in range(mpar): 203 if index2*mpar+j >= index1*mpar+i: 204 cov_subset_vals.append(cov[index1*mpar+i, 205 index2*mpar+j]) 206 207 cov_subset = np.zeros([sub_cov_dim, sub_cov_dim]) 208 tri_indices = np.triu_indices(sub_cov_dim) 209 cov_subset[tri_indices] = cov_subset_vals 210 211 res_cov_range.file2res.parameters = parameters[mask] 212 res_cov_range.covariance = cov_subset 213 return res_cov_range 214 215 def sample(self, n_samples): 216 """Sample resonance parameters based on the covariances provided 217 within an ENDF evaluation. 218 219 Parameters 220 ---------- 221 n_samples : int 222 The number of samples to produce 223 224 Returns 225 ------- 226 samples : list of openmc.data.ResonanceCovarianceRange objects 227 List of samples size `n_samples` 228 229 """ 230 warn_str = 'Sampling routine does not guarantee positive values for '\ 231 'parameters. This can lead to undefined behavior in the '\ 232 'reconstruction routine.' 233 warnings.warn(warn_str) 234 parameters = self.parameters 235 cov = self.covariance 236 237 # Symmetrizing covariance matrix 238 cov = cov + cov.T - np.diag(cov.diagonal()) 239 formalism = self.formalism 240 mpar = self.mpar 241 samples = [] 242 243 # Handling MLBW/SLBW sampling 244 if formalism == 'mlbw' or formalism == 'slbw': 245 params = ['energy', 'neutronWidth', 'captureWidth', 'fissionWidth', 246 'competitiveWidth'] 247 param_list = params[:mpar] 248 mean_array = parameters[param_list].values 249 mean = mean_array.flatten() 250 par_samples = np.random.multivariate_normal(mean, cov, 251 size=n_samples) 252 spin = parameters['J'].values 253 l_value = parameters['L'].values 254 for sample in par_samples: 255 energy = sample[0::mpar] 256 gn = sample[1::mpar] 257 gg = sample[2::mpar] 258 gf = sample[3::mpar] if mpar > 3 else parameters['fissionWidth'].values 259 gx = sample[4::mpar] if mpar > 4 else parameters['competitiveWidth'].values 260 gt = gn + gg + gf + gx 261 262 records = [] 263 for j, E in enumerate(energy): 264 records.append([energy[j], l_value[j], spin[j], gt[j], 265 gn[j], gg[j], gf[j], gx[j]]) 266 columns = ['energy', 'L', 'J', 'totalWidth', 'neutronWidth', 267 'captureWidth', 'fissionWidth', 'competitiveWidth'] 268 sample_params = pd.DataFrame.from_records(records, 269 columns=columns) 270 # Copy ResonanceRange object 271 res_range = copy.copy(self.file2res) 272 res_range.parameters = sample_params 273 samples.append(res_range) 274 275 # Handling RM sampling 276 elif formalism == 'rm': 277 params = ['energy', 'neutronWidth', 'captureWidth', 278 'fissionWidthA', 'fissionWidthB'] 279 param_list = params[:mpar] 280 mean_array = parameters[param_list].values 281 mean = mean_array.flatten() 282 par_samples = np.random.multivariate_normal(mean, cov, 283 size=n_samples) 284 spin = parameters['J'].values 285 l_value = parameters['L'].values 286 for sample in par_samples: 287 energy = sample[0::mpar] 288 gn = sample[1::mpar] 289 gg = sample[2::mpar] 290 gfa = sample[3::mpar] if mpar > 3 else parameters['fissionWidthA'].values 291 gfb = sample[4::mpar] if mpar > 3 else parameters['fissionWidthB'].values 292 293 records = [] 294 for j, E in enumerate(energy): 295 records.append([energy[j], l_value[j], spin[j], gn[j], 296 gg[j], gfa[j], gfb[j]]) 297 columns = ['energy', 'L', 'J', 'neutronWidth', 298 'captureWidth', 'fissionWidthA', 'fissionWidthB'] 299 sample_params = pd.DataFrame.from_records(records, 300 columns=columns) 301 # Copy ResonanceRange object 302 res_range = copy.copy(self.file2res) 303 res_range.parameters = sample_params 304 samples.append(res_range) 305 306 return samples 307 308 309class MultiLevelBreitWignerCovariance(ResonanceCovarianceRange): 310 """Multi-level Breit-Wigner resolved resonance formalism covariance data. 311 Parameters 312 ---------- 313 energy_min : float 314 Minimum energy of the resolved resonance range in eV 315 energy_max : float 316 Maximum energy of the resolved resonance range in eV 317 318 Attributes 319 ---------- 320 energy_min : float 321 Minimum energy of the resolved resonance range in eV 322 energy_max : float 323 Maximum energy of the resolved resonance range in eV 324 parameters : pandas.DataFrame 325 Resonance parameters 326 covariance : numpy.array 327 The covariance matrix contained within the ENDF evaluation 328 mpar : int 329 Number of parameters in covariance matrix for each individual resonance 330 lcomp : int 331 Flag indicating format of the covariance matrix within the ENDF file 332 file2res : openmc.data.ResonanceRange object 333 Corresponding resonance range with File 2 data. 334 formalism : str 335 String descriptor of formalism 336 337 """ 338 339 def __init__(self, energy_min, energy_max, parameters, covariance, mpar, 340 lcomp, file2res): 341 super().__init__(energy_min, energy_max) 342 self.parameters = parameters 343 self.covariance = covariance 344 self.mpar = mpar 345 self.lcomp = lcomp 346 self.file2res = copy.copy(file2res) 347 self.formalism = 'mlbw' 348 349 @classmethod 350 def from_endf(cls, ev, file_obj, items, resonance): 351 """Create MLBW covariance data from an ENDF evaluation. 352 353 Parameters 354 ---------- 355 ev : openmc.data.endf.Evaluation 356 ENDF evaluation 357 file_obj : file-like object 358 ENDF file positioned at the second record of a resonance range 359 subsection in MF=32, MT=151 360 items : list 361 Items from the CONT record at the start of the resonance range 362 subsection 363 resonance : openmc.data.ResonanceRange object 364 Corresponding resonance range with File 2 data. 365 366 Returns 367 ------- 368 openmc.data.MultiLevelBreitWignerCovariance 369 Multi-level Breit-Wigner resonance covariance parameters 370 371 """ 372 373 # Read energy-dependent scattering radius if present 374 energy_min, energy_max = items[0:2] 375 nro, naps = items[4:6] 376 if nro != 0: 377 params, ape = endf.get_tab1_record(file_obj) 378 379 # Other scatter radius parameters 380 items = endf.get_cont_record(file_obj) 381 target_spin = items[0] 382 lcomp = items[3] # Flag for compatibility 0, 1, 2 - 2 is compact form 383 nls = items[4] # number of l-values 384 385 # Build covariance matrix for General Resolved Resonance Formats 386 if lcomp == 1: 387 items = endf.get_cont_record(file_obj) 388 # Number of short range type resonance covariances 389 num_short_range = items[4] 390 # Number of long range type resonance covariances 391 num_long_range = items[5] 392 393 # Read resonance widths, J values, etc 394 records = [] 395 for i in range(num_short_range): 396 items, values = endf.get_list_record(file_obj) 397 mpar = items[2] 398 num_res = items[5] 399 num_par_vals = num_res*6 400 res_values = values[:num_par_vals] 401 cov_values = values[num_par_vals:] 402 403 energy = res_values[0::6] 404 spin = res_values[1::6] 405 gt = res_values[2::6] 406 gn = res_values[3::6] 407 gg = res_values[4::6] 408 gf = res_values[5::6] 409 410 for i, E in enumerate(energy): 411 records.append([energy[i], spin[i], gt[i], gn[i], 412 gg[i], gf[i]]) 413 414 # Build the upper-triangular covariance matrix 415 cov_dim = mpar*num_res 416 cov = np.zeros([cov_dim, cov_dim]) 417 indices = np.triu_indices(cov_dim) 418 cov[indices] = cov_values 419 420 # Compact format - Resonances and individual uncertainties followed by 421 # compact correlations 422 elif lcomp == 2: 423 items, values = endf.get_list_record(file_obj) 424 mean = items 425 num_res = items[5] 426 energy = values[0::12] 427 spin = values[1::12] 428 gt = values[2::12] 429 gn = values[3::12] 430 gg = values[4::12] 431 gf = values[5::12] 432 par_unc = [] 433 for i in range(num_res): 434 res_unc = values[i*12+6 : i*12+12] 435 # Delete 0 values (not provided, no fission width) 436 # DAJ/DGT always zero, DGF sometimes nonzero [1, 2, 5] 437 res_unc_nonzero = [] 438 for j in range(6): 439 if j in [1, 2, 5] and res_unc[j] != 0.0: 440 res_unc_nonzero.append(res_unc[j]) 441 elif j in [0, 3, 4]: 442 res_unc_nonzero.append(res_unc[j]) 443 par_unc.extend(res_unc_nonzero) 444 445 records = [] 446 for i, E in enumerate(energy): 447 records.append([energy[i], spin[i], gt[i], gn[i], 448 gg[i], gf[i]]) 449 450 corr = endf.get_intg_record(file_obj) 451 cov = np.diag(par_unc).dot(corr).dot(np.diag(par_unc)) 452 453 # Compatible resolved resonance format 454 elif lcomp == 0: 455 cov = np.zeros([4, 4]) 456 records = [] 457 cov_index = 0 458 for i in range(nls): 459 items, values = endf.get_list_record(file_obj) 460 num_res = items[5] 461 for j in range(num_res): 462 one_res = values[18*j:18*(j+1)] 463 res_values = one_res[:6] 464 cov_values = one_res[6:] 465 records.append(list(res_values)) 466 467 # Populate the coviariance matrix for this resonance 468 # There are no covariances between resonances in lcomp=0 469 cov[cov_index, cov_index] = cov_values[0] 470 cov[cov_index+1, cov_index+1 : cov_index+2] = cov_values[1:2] 471 cov[cov_index+1, cov_index+3] = cov_values[4] 472 cov[cov_index+2, cov_index+2] = cov_values[3] 473 cov[cov_index+2, cov_index+3] = cov_values[5] 474 cov[cov_index+3, cov_index+3] = cov_values[6] 475 476 cov_index += 4 477 if j < num_res-1: # Pad matrix for additional values 478 cov = np.pad(cov, ((0, 4), (0, 4)), 'constant', 479 constant_values=0) 480 481 # Create pandas DataFrame with resonance data, currently 482 # redundant with data.IncidentNeutron.resonance 483 columns = ['energy', 'J', 'totalWidth', 'neutronWidth', 484 'captureWidth', 'fissionWidth'] 485 parameters = pd.DataFrame.from_records(records, columns=columns) 486 # Determine mpar (number of parameters for each resonance in 487 # covariance matrix) 488 nparams, params = parameters.shape 489 covsize = cov.shape[0] 490 mpar = int(covsize/nparams) 491 # Add parameters from File 2 492 parameters = _add_file2_contributions(parameters, 493 resonance.parameters) 494 # Create instance of class 495 mlbw = cls(energy_min, energy_max, parameters, cov, mpar, lcomp, 496 resonance) 497 return mlbw 498 499 500class SingleLevelBreitWignerCovariance(MultiLevelBreitWignerCovariance): 501 """Single-level Breit-Wigner resolved resonance formalism covariance data. 502 Single-level Breit-Wigner resolved resonance data is is identified by LRF=1 503 in the ENDF-6 format. 504 505 Parameters 506 ---------- 507 energy_min : float 508 Minimum energy of the resolved resonance range in eV 509 energy_max : float 510 Maximum energy of the resolved resonance range in eV 511 512 Attributes 513 ---------- 514 energy_min : float 515 Minimum energy of the resolved resonance range in eV 516 energy_max : float 517 Maximum energy of the resolved resonance range in eV 518 parameters : pandas.DataFrame 519 Resonance parameters 520 covariance : numpy.array 521 The covariance matrix contained within the ENDF evaluation 522 mpar : int 523 Number of parameters in covariance matrix for each individual resonance 524 formalism : str 525 String descriptor of formalism 526 lcomp : int 527 Flag indicating format of the covariance matrix within the ENDF file 528 file2res : openmc.data.ResonanceRange object 529 Corresponding resonance range with File 2 data. 530 """ 531 532 def __init__(self, energy_min, energy_max, parameters, covariance, mpar, 533 lcomp, file2res): 534 super().__init__(energy_min, energy_max, parameters, covariance, mpar, 535 lcomp, file2res) 536 self.formalism = 'slbw' 537 538 539class ReichMooreCovariance(ResonanceCovarianceRange): 540 """Reich-Moore resolved resonance formalism covariance data. 541 542 Reich-Moore resolved resonance data is identified by LRF=3 in the ENDF-6 543 format. 544 545 Parameters 546 ---------- 547 energy_min : float 548 Minimum energy of the resolved resonance range in eV 549 energy_max : float 550 Maximum energy of the resolved resonance range in eV 551 552 Attributes 553 ---------- 554 energy_min : float 555 Minimum energy of the resolved resonance range in eV 556 energy_max : float 557 Maximum energy of the resolved resonance range in eV 558 parameters : pandas.DataFrame 559 Resonance parameters 560 covariance : numpy.array 561 The covariance matrix contained within the ENDF evaluation 562 lcomp : int 563 Flag indicating format of the covariance matrix within the ENDF file 564 mpar : int 565 Number of parameters in covariance matrix for each individual resonance 566 file2res : openmc.data.ResonanceRange object 567 Corresponding resonance range with File 2 data. 568 formalism : str 569 String descriptor of formalism 570 """ 571 572 def __init__(self, energy_min, energy_max, parameters, covariance, mpar, 573 lcomp, file2res): 574 super().__init__(energy_min, energy_max) 575 self.parameters = parameters 576 self.covariance = covariance 577 self.mpar = mpar 578 self.lcomp = lcomp 579 self.file2res = copy.copy(file2res) 580 self.formalism = 'rm' 581 582 @classmethod 583 def from_endf(cls, ev, file_obj, items, resonance): 584 """Create Reich-Moore resonance covariance data from an ENDF 585 evaluation. Includes the resonance parameters contained separately in 586 File 32. 587 588 Parameters 589 ---------- 590 ev : openmc.data.endf.Evaluation 591 ENDF evaluation 592 file_obj : file-like object 593 ENDF file positioned at the second record of a resonance range 594 subsection in MF=2, MT=151 595 items : list 596 Items from the CONT record at the start of the resonance range 597 subsection 598 resonance : openmc.data.Resonance object 599 openmc.data.Resonanance object generated from the same evaluation 600 used to import values not contained in File 32 601 602 Returns 603 ------- 604 openmc.data.ReichMooreCovariance 605 Reich-Moore resonance covariance parameters 606 607 """ 608 # Read energy-dependent scattering radius if present 609 energy_min, energy_max = items[0:2] 610 nro, naps = items[4:6] 611 if nro != 0: 612 params, ape = endf.get_tab1_record(file_obj) 613 614 # Other scatter radius parameters 615 items = endf.get_cont_record(file_obj) 616 target_spin = items[0] 617 lcomp = items[3] # Flag for compatibility 0, 1, 2 - 2 is compact form 618 nls = items[4] # Number of l-values 619 620 # Build covariance matrix for General Resolved Resonance Formats 621 if lcomp == 1: 622 items = endf.get_cont_record(file_obj) 623 # Number of short range type resonance covariances 624 num_short_range = items[4] 625 # Number of long range type resonance covariances 626 num_long_range = items[5] 627 # Read resonance widths, J values, etc 628 channel_radius = {} 629 scattering_radius = {} 630 records = [] 631 for i in range(num_short_range): 632 items, values = endf.get_list_record(file_obj) 633 mpar = items[2] 634 num_res = items[5] 635 num_par_vals = num_res*6 636 res_values = values[:num_par_vals] 637 cov_values = values[num_par_vals:] 638 639 energy = res_values[0::6] 640 spin = res_values[1::6] 641 gn = res_values[2::6] 642 gg = res_values[3::6] 643 gfa = res_values[4::6] 644 gfb = res_values[5::6] 645 646 for i, E in enumerate(energy): 647 records.append([energy[i], spin[i], gn[i], gg[i], 648 gfa[i], gfb[i]]) 649 650 # Build the upper-triangular covariance matrix 651 cov_dim = mpar*num_res 652 cov = np.zeros([cov_dim, cov_dim]) 653 indices = np.triu_indices(cov_dim) 654 cov[indices] = cov_values 655 656 # Compact format - Resonances and individual uncertainties followed by 657 # compact correlations 658 elif lcomp == 2: 659 items, values = endf.get_list_record(file_obj) 660 num_res = items[5] 661 energy = values[0::12] 662 spin = values[1::12] 663 gn = values[2::12] 664 gg = values[3::12] 665 gfa = values[4::12] 666 gfb = values[5::12] 667 par_unc = [] 668 for i in range(num_res): 669 res_unc = values[i*12+6 : i*12+12] 670 # Delete 0 values (not provided in evaluation) 671 res_unc = [x for x in res_unc if x != 0.0] 672 par_unc.extend(res_unc) 673 674 records = [] 675 for i, E in enumerate(energy): 676 records.append([energy[i], spin[i], gn[i], gg[i], 677 gfa[i], gfb[i]]) 678 679 corr = endf.get_intg_record(file_obj) 680 cov = np.diag(par_unc).dot(corr).dot(np.diag(par_unc)) 681 682 # Create pandas DataFrame with resonacne data 683 columns = ['energy', 'J', 'neutronWidth', 'captureWidth', 684 'fissionWidthA', 'fissionWidthB'] 685 parameters = pd.DataFrame.from_records(records, columns=columns) 686 687 # Determine mpar (number of parameters for each resonance in 688 # covariance matrix) 689 nparams, params = parameters.shape 690 covsize = cov.shape[0] 691 mpar = int(covsize/nparams) 692 693 # Add parameters from File 2 694 parameters = _add_file2_contributions(parameters, 695 resonance.parameters) 696 # Create instance of ReichMooreCovariance 697 rmc = cls(energy_min, energy_max, parameters, cov, mpar, lcomp, 698 resonance) 699 return rmc 700 701 702_FORMALISMS = { 703 0: ResonanceCovarianceRange, 704 1: SingleLevelBreitWignerCovariance, 705 2: MultiLevelBreitWignerCovariance, 706 3: ReichMooreCovariance 707 # 7: RMatrixLimitedCovariance 708} 709