Implement vector fitting to replace external vectfit package

[azimgivron]

Description

This pull request removes the dependency on the unmaintained external vectfit package (https://github.com/liangjg/vectfit) used in OpenMC's windowed multipole (WMP) infrastructure. In its place, a vector fitting is introduced based on the vectorFitting.py module from scikit-rf project.

The implementation will be made in several steps as suggested by Paul:

• Confirm whether the wrapper code it produced is correct or needs updating
• Vendor a stripped version of the necessary vector fitting functionality from scikit-rf's vectorFitting.py module
• Test it out by generating multipole data from an ENDF file and compare what you get with the updated implementation compared to the vectfit module. Hopefully it agrees but if not it will require some more digging.

This change improves long-term maintainability and compatibility of OpenMC.

Fixes #3487

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 15) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[azimgivron]

This PR is currently a draft. The code is a copy-pasting of the suggestion by chatGPT provided by Paul.
The following repository provides an installation of openmc with vectfit dependency as it currently is implemented:
https://github.com/azimgivron/openmc-vectfit.git
This will be used to make the comparisons between the implementations.

[azimgivron]

Hello guys,

I've taken the time to re-implement the code using NumPy and SciPy. The unit tests are adapted from Jingang Liang’s implementation for consistency.

To validate the results, I generated WMP data from ENDF/B-VII for both Na-23 and Fe-56, and ran both implementations on these datasets. I will soon share plots comparing the absolute error in absorption cross sections between the two codes but there almost 0.

I also benchmarked the compute time, though I haven’t yet post-processed those results. In any case, I believe this additional performance data isn’t critical for this PR.

[azimgivron]

📊 Interpolation Benchmark: C++ vs Python

Note that the legend says $|\text{C++}-\text{Python}|$ where it should say $\frac{\text{C++}-\text{Python}}{\text{C++}}$.

This benchmark compares the C++ and Python implementations of vectfit for the WMP absorption cross sections of Na23 and Fe56. The relative error plots show that the results are close but differ at some points, I guess its energy frontier between interpolations, though further investigation could be done.

As expected, the Python implementation is slightly slower than the C++ counterpart, though still functionally equivalent in output.

🧬 Interpolation Metadata

Na23:
- 56 resonance poles
- 425 energy windows
- 0.63 poles per window
- 5th-order Chebyshev fit for the smooth background

Fe56:
- 408 resonance poles
- 1273 energy windows
- 2.7 poles per window
- 2nd-order Chebyshev fit for the smooth background

[paulromano]

@azimgivron Thanks so much for going through this conversion! I made a few small updates on your branch but otherwise it looks pretty good. I do have one important question about a scaling below (wondering if this is affecting the results) below:

[azimgivron]

Hi, is this PR blocked because something is expected from me ?

[paulromano]

@azimgivron Any chance you can update the results you had shared above with the normalization correction? I'm curious how much of the relative error shown before may be due to that.

[GuySten]

This PR looks ready to be merge for me.
@azimgivron, could you provide the scripts used to make the comparison between python and cpp vector fitting?

[GuySten]

Looks good to me.

[GuySten]

I also compared the cpp and python wmp absorption cross sections:

The cpp version has:
Na23:

• 58 resonance poles
• 320 energy windows
• 0.784375 poles per window
• 5th-order Chebyshev fit for the smooth background
  Fe56:
• 407 resonance poles
• 1273 energy windows
• 2.633935585231736 poles per window
• 2th-order Chebyshev fit for the smooth background

while the python version has:
Na23:

• 60 resonance poles
• 320 energy windows
• 0.6 poles per window
• 5th-order Chebyshev fit for the smooth background
  Fe56:
• 403 resonance poles
• 1273 energy windows
• 2.6441476826394346 poles per window
• 2th-order Chebyshev fit for the smooth background.

The comparison is pretty good for Na23 (under 1%).
But is a bit worse for Fe56 (around 3%).
We might be able to improve that if we want by tightening the convergence tolerances of the vector fitting method.

If no one objects I am planning to merge this feature soon.

[paulromano]

A relative error of 3% in the Fe56 cross section seems unacceptable to me, especially with no specific explanation of why they are different. Another way of looking at this would be to ask how close the WMP-evaluated cross section is relative to the original tabular data that it is fitting. As far as I'm aware, it should do much better than 3%. Looking at this paper, vector fitting should be able to match cross sections within 0.1%.

[GuySten]

I further compared the python and cpp to the original endf evaluation and got:

So the error while using this feature got smaller.

[GuySten]

@paulromano, what do you think should be done?

According to my testing the vector fitting that we already have more than 5% error and by using the python version it get reduced to around 3%.
Should we merge it as is and open an issue that the default vector fitting error is too large?

I think that the fact that the vector fitting is done currently using cpp code that without the docker image @azimgivron made I could not manage to install is a big barrier for users. So I would not be surprised that it has bugs no one noticed.

EDIT:
You can also see that both the python and cpp implemantation has cross sections that do not reach 20MeV.
So the WMP feature is still not mature enough.

[paulromano]

That's very interesting. I didn't realize that the original was that far off. @bforget do you have any ideas on why the original implementation that @liangjg put together would be off by that much?

Regarding cross sections not going up to 20 MeV -- I believe that is by design as the vector fitting is only needed over the resolved resonance range. At fast energies it is assumed there is effectively no temperature dependence (which obviates the need for WMP data).

In any case, given the data you've shared, I would be more amenable to getting this merged and tagging this as an issue to be looked into further as to why the comparison with ENDF data doesn't look great. I'll give the PR a final look as soon as I can.

[paulromano]

I made a few small changes here and I think this is good to go now. I'll leave this open until tomorrow in case there are any final comments. Many thanks to you both @azimgivron and @GuySten!