diff --git a/CHANGELOG.md b/CHANGELOG.md index 38ac9f48..7f4a56f6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,7 +9,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added -- Add an interface to the Wolfram Language to allow Rust APIs to be called in +- Added proper LHAPDF drop-in compatibility layer for no-code migration. +- Added an interface to the Wolfram Language to allow Rust APIs to be called in Mathematica. - Added the logic to determine the Euclidean distance of a point to the closest subgrid in order to allow extrapolation. diff --git a/neopdf_capi/src/include/NeoPDF.hpp b/neopdf_capi/src/include/NeoPDF.hpp index 66d0c79b..a3de2bc9 100644 --- a/neopdf_capi/src/include/NeoPDF.hpp +++ b/neopdf_capi/src/include/NeoPDF.hpp @@ -454,3 +454,44 @@ class GridWriter { }; } // namespace neopdf + +/** @brief LHAPDF compatibility for no-code migration. */ +namespace NEOLHAPDF { + +class PDF { + protected: + std::unique_ptr _neopdf; + + PDF() = default; + explicit PDF(std::unique_ptr&& pdf) : _neopdf(std::move(pdf)) {} + + public: + virtual ~PDF() = default; + + virtual double xfxQ2(int id, double x, double q2) const { + return _neopdf->xfxQ2(id, x, q2); + } + + virtual double alphasQ2(double q2) const { + return _neopdf->alphasQ2(q2); + } + + double xMin() const { return _neopdf->x_min(); } + double xMax() const { return _neopdf->x_max(); } + double q2Min() const { return _neopdf->q2_min(); } + double q2Max() const { return _neopdf->q2_max(); } +}; + +class GridPDF : public PDF { + public: + explicit GridPDF(std::unique_ptr&& pdf) : PDF(std::move(pdf)) {} +}; + +inline PDF* mkPDF(const std::string& name, int member = 0) { + std::unique_ptr neopdf_ptr(new neopdf::NeoPDF(name, member)); + return new GridPDF(std::move(neopdf_ptr)); +} + +inline void setVerbosity(int /*verbosity*/) { } + +} // namespace LHAPDF diff --git a/neopdf_capi/src/lib.rs b/neopdf_capi/src/lib.rs index d4b590a9..a50d73ea 100644 --- a/neopdf_capi/src/lib.rs +++ b/neopdf_capi/src/lib.rs @@ -10,6 +10,8 @@ use neopdf::parser::SubgridData; use neopdf::pdf::PDF; use neopdf::writer::GridArrayCollection; +const DEFAULT_PIDS: [i32; 14] = [21, -6, -5, -4, -3, -2, -1, 1, 2, 3, 4, 5, 6, 22]; + /// Result codes for `NeoPDF` operations #[repr(C)] #[derive(Debug, Clone, Copy, PartialEq, Eq)] @@ -1073,3 +1075,237 @@ pub unsafe extern "C" fn neopdf_grid_compress( Err(_) => NeopdfResult::ErrorMemoryError, } } + +// LHAPDF C-API drop-in compatibility layer. +/// +/// This global state stores the loaded PDF set and the currently selected +/// member index. It is not thread-safe and mirrors the behavior of legacy +/// LHAPDF interfaces. +struct LhapdfState { + pdf_set: Option>, + member: usize, +} + +static mut LHAPDF_STATE: LhapdfState = LhapdfState { + pdf_set: None, + member: 0, +}; + +/// Sets LHAPDF runtime parameters from a string (no-op). +/// +/// This function is provided for compatibility with the LHAPDF C API and does +/// nothing in this implementation. +/// +/// # Safety +/// +/// The pointer `line` must be a valid, null-terminated C string if non-null. +#[no_mangle] +pub const unsafe extern "C" fn setlhaparm(_line: *const c_char) {} + +/// Fortran name-mangled variant of `setlhaparm` (no-op). +/// +/// This function is provided for compatibility with the LHAPDF Fortran API and +/// does nothing in this implementation. +/// +/// # Safety +/// +/// The pointer `line` must be valid for reads of `len` bytes if non-null. +#[no_mangle] +pub const unsafe extern "C" fn setlhaparm_(_line: *const c_char, _len: isize) {} + +/// Initializes a PDF set by its name/path and loads all members. +/// +/// The loaded set is stored in a global state used by the other LHAPDF-compatible +/// functions, with the current member index reset to 0. +/// +/// # Panics +/// +/// Panics if `name` is not valid UTF-8 or is not a valid, null-terminated C string. +/// +/// # Safety +/// +/// `name` must be a valid, null-terminated C string pointing to a UTF-8 path or set name. +#[no_mangle] +pub unsafe extern "C" fn initpdfsetbyname(name: *const c_char) { + unsafe { + let c_str = CStr::from_ptr(name); + let pdf_name = c_str.to_str().expect("Invalid UTF-8 string"); + let pdfs = PDF::load_pdfs(pdf_name); + LHAPDF_STATE.pdf_set = Some(pdfs); + LHAPDF_STATE.member = 0; + } +} + +/// Fortran name-mangled variant of `initpdfsetbyname`. +/// +/// Reads a fixed-length Fortran character buffer, trims trailing spaces, and loads +/// the corresponding PDF set into the global state. The current member index is +/// reset to 0. +/// +/// # Panics +/// +/// TODO +/// +/// # Safety +/// +/// `name` must be valid for reads of `len` bytes. The buffer may not be +/// null-terminated; trailing spaces are trimmed. +#[no_mangle] +#[allow(clippy::cast_sign_loss)] +pub unsafe extern "C" fn initpdfsetbyname_(name: *const c_char, len: c_int) { + unsafe { + let name_slice = slice::from_raw_parts(name.cast::(), len as usize); + let pdf_name = std::str::from_utf8(name_slice).unwrap().trim_end(); + let pdfs = PDF::load_pdfs(pdf_name); + LHAPDF_STATE.pdf_set = Some(pdfs); + LHAPDF_STATE.member = 0; + } +} + +/// Selects the active member of the currently loaded PDF set. +/// +/// # Safety +/// +/// This function does not perform bounds checking here; subsequent calls that +/// use the member will return early if the index is out of range. +#[no_mangle] +#[allow(clippy::cast_sign_loss)] +pub unsafe extern "C" fn initpdf(member: c_int) { + unsafe { LHAPDF_STATE.member = member as usize }; +} + +/// Fortran name-mangled variant of `initpdf`. +/// +/// # Safety +/// +/// `member` must be a valid pointer to an integer. +#[no_mangle] +#[allow(clippy::cast_sign_loss)] +pub unsafe extern "C" fn initpdf_(member: *const c_int) { + unsafe { LHAPDF_STATE.member = *member as usize }; +} + +/// Evaluates parton distribution functions at given `(x, q)` for the active member. +/// +/// # Safety +/// +/// - `f` must point to writable memory for at least 13 `c_double` values. +/// - Requires that a PDF set has been initialized via `initpdfsetbyname` or its +/// Fortran variant. If no set is loaded or the member index is out of range, +/// the function returns without writing. +#[no_mangle] +pub unsafe extern "C" fn evolvepdf(x: c_double, q: c_double, f: *mut c_double) { + unsafe { + let state_ptr = &raw const LHAPDF_STATE; + let pdf_set_ptr = &raw const (*state_ptr).pdf_set; + + if let Some(pdfs) = &*pdf_set_ptr { + let member = (*state_ptr).member; + if member >= pdfs.len() { + return; + } + let pdf = &pdfs[member]; + let q2 = q * q; + + let (available_pids, _) = pdf.pids().clone().into_raw_vec_and_offset(); + let out_slice = slice::from_raw_parts_mut(f, 14); + + for (out, &pid) in out_slice.iter_mut().zip(DEFAULT_PIDS.iter()) { + *out = if available_pids.contains(&pid) { + pdf.xfxq2(pid, &[x, q2]) + } else { + 0.0 + }; + } + } + } +} + +/// Fortran name-mangled variant of `evolvepdf`. +/// +/// # Safety +/// +/// - `x`, `q` must be valid pointers to `c_double` values. +/// - `f` must point to writable memory for at least 13 `c_double` values. +/// - A PDF set must have been initialized and the member index must be valid or +/// the function will return without writing. +#[no_mangle] +pub unsafe extern "C" fn evolvepdf_(x: *const c_double, q: *const c_double, f: *mut c_double) { + unsafe { + let pdf_set_ptr = &raw const LHAPDF_STATE.pdf_set; + + if let Some(pdfs) = &*pdf_set_ptr { + let member = (&raw const LHAPDF_STATE.member).read(); + if member >= pdfs.len() { + return; + } + let pdf = &pdfs[member]; + let q2 = (*q) * (*q); + + let (available_pids, _) = pdf.pids().clone().into_raw_vec_and_offset(); + let out_slice = slice::from_raw_parts_mut(f, 14); + + for (out, &pid) in out_slice.iter_mut().zip(DEFAULT_PIDS.iter()) { + *out = if available_pids.contains(&pid) { + pdf.xfxq2(pid, &[*x, q2]) + } else { + 0.0 + }; + } + } + } +} + +/// Evaluates the strong coupling `alpha_s` at scale `q` for the active member. +/// +/// # Safety +/// +/// - Requires that a PDF set has been initialized via `initpdfsetbyname` or its +/// Fortran variant. If no set is loaded or the member index is out of range, +/// the function returns 0.0. +#[no_mangle] +pub unsafe extern "C" fn alphaspdf(q: c_double) -> c_double { + unsafe { + let state_ptr = &raw const LHAPDF_STATE; + let pdf_set_ptr = &raw const (*state_ptr).pdf_set; + + if let Some(pdfs) = &*pdf_set_ptr { + let member = (*state_ptr).member; + if member >= pdfs.len() { + return 0.0; + } + let pdf = &pdfs[member]; + let q2 = q * q; + pdf.alphas_q2(q2) + } else { + 0.0 + } + } +} + +/// Fortran name-mangled variant of `alphaspdf`. +/// +/// # Safety +/// +/// - `q` must be a valid pointer to a `c_double` value. +/// - A PDF set must have been initialized and the member index must be valid or +/// the function will return 0.0. +#[no_mangle] +pub unsafe extern "C" fn alphaspdf_(q: *const c_double) -> c_double { + unsafe { + let state_ptr = &raw const LHAPDF_STATE; + let pdf_set_ptr = &raw const (*state_ptr).pdf_set; + + if let Some(pdfs) = &*pdf_set_ptr { + let member = (*state_ptr).member; + if member >= pdfs.len() { + return 0.0; + } + let pdf = &pdfs[member]; + let q2 = (*q) * (*q); + pdf.alphas_q2(q2) + } else { + 0.0 + } + } +} diff --git a/neopdf_capi/tests/Makefile b/neopdf_capi/tests/Makefile index 4024f217..46eb724f 100644 --- a/neopdf_capi/tests/Makefile +++ b/neopdf_capi/tests/Makefile @@ -5,7 +5,7 @@ NEOPDF_DEPS != pkg-config --cflags --libs neopdf_capi LHAPDF_DEPS != pkg-config --cflags --libs lhapdf MATH_LIBS = -lm -PROGRAMS = check-capi check-oop check-writer check-writer-oop check-xapi check-xwriter +PROGRAMS = check-capi check-oop check-writer check-writer-oop check-xapi check-xwriter check-lhapdf-compatibility all: $(PROGRAMS) @@ -30,6 +30,9 @@ check-xapi: check-capi.c check-xwriter: check-writer.c $(CC) $(CXXFLAGS) $< $(NEOPDF_DEPS) $(MATH_LIBS) -o $@ +check-lhapdf-compatibility: check-lhapdf-compatibility.cpp + $(CXX) $(CXXFLAGS) $< $(NEOPDF_DEPS) -o $@ + .PHONY: clean clean: diff --git a/neopdf_capi/tests/check-lhapdf-compatibility.cpp b/neopdf_capi/tests/check-lhapdf-compatibility.cpp new file mode 100644 index 00000000..496c3c55 --- /dev/null +++ b/neopdf_capi/tests/check-lhapdf-compatibility.cpp @@ -0,0 +1,99 @@ +#include +#include +#include +#include +#include +#include + +using namespace NEOLHAPDF; + +const double TOLERANCE = 1e-12; + +void test_lhapdf_compatibility_oop() { + std::cout << "=== Test LHAPDF C++ OOP Compatibility Layer ===\n"; + + std::string pdfname = "NNPDF40_nnlo_as_01180"; + + // Use the compatibility layer + setVerbosity(0); + PDF* neo_lha_pdf = mkPDF(pdfname, 0); + + // Use the native NeoPDF C API for comparison + NeoPDFWrapper* neo_pdf = neopdf_pdf_load(pdfname.c_str(), 0); + + double x = 1e-4; + double q2 = 100.0; + int pid = 21; + + double result = neo_lha_pdf->xfxQ2(pid, x, q2); + double expected = neopdf_pdf_xfxq2(neo_pdf, pid, x, q2); + double reldif = std::abs(result - expected) / expected; + + std::cout << "pid=" << pid << ", x=" << x << ", Q2=" << q2 << std::endl; + std::cout << "Native NeoPDF: " << expected << std::endl; + std::cout << "NeoPDF (LHAPDF compat): " << result << std::endl; + std::cout << "Relative difference: " << reldif << std::endl; + + double as_result = neo_lha_pdf->alphasQ2(q2); + double as_expected = neopdf_pdf_alphas_q2(neo_pdf, q2); + double as_reldif = std::abs(as_result - as_expected) / as_expected; + + std::cout << "Q2=" << q2 << std::endl; + std::cout << "Native NeoPDF: " << as_expected << std::endl; + std::cout << "NeoPDF (LHAPDF compat): " << as_result << std::endl; + std::cout << "Relative difference: " << as_reldif << std::endl; + + assert(reldif < TOLERANCE); + assert(as_reldif < TOLERANCE); + + delete neo_lha_pdf; + neopdf_pdf_free(neo_pdf); + + std::cout << "LHAPDF C++ OOP Compatibility test passed." << std::endl; +} + +void test_lhapdf_compatibility_c() { + printf("=== Test LHAPDF C Compatibility Layer ===\n"); + + const char* pdfname = "NNPDF40_nnlo_as_01180"; + + // Use the compatibility layer + initpdfsetbyname(pdfname); + initpdf(0); + double x = 0.1; + double q = 10.0; + double xfxs[14]; + evolvepdf(x, q, xfxs); + double a_s = alphaspdf(q); + + // Use the native NeoPDF API for comparison + NeoPDFWrapper* neo_pdf = neopdf_pdf_load(pdfname, 0); + double expected_g = neopdf_pdf_xfxq2(neo_pdf, -1, x, q*q); + neopdf_pdf_free(neo_pdf); + double as_expected = neopdf_pdf_alphas_q2(neo_pdf, q*q); + + double reldiff = fabs(xfxs[6] - expected_g) / expected_g; + double as_reldiff = fabs(a_s - as_expected) / as_expected; + + printf("Comparing gluon PDF at x=0.1, Q=10.0\n"); + printf("Native NeoPDF: %e\n", expected_g); + printf("NeoPDF (LHAPDF compat): %e\n", xfxs[6]); + printf("Relative difference: %e\n", reldiff); + + printf("Comparing AlphaS at Q=10.0\n"); + printf("Native NeoPDF: %e\n", as_expected); + printf("NeoPDF (LHAPDF compat): %e\n", a_s); + printf("Relative difference: %e\n", as_reldiff); + + assert(reldiff < TOLERANCE); + assert(as_reldiff < TOLERANCE); + + printf("LHAPDF C Compatibility test passed.\n"); +} + +int main() { + test_lhapdf_compatibility_oop(); + test_lhapdf_compatibility_c(); + + return 0; +} diff --git a/neopdf_capi/tests/check-lhapdf-compatibility.output b/neopdf_capi/tests/check-lhapdf-compatibility.output new file mode 100644 index 00000000..7b8c4a53 --- /dev/null +++ b/neopdf_capi/tests/check-lhapdf-compatibility.output @@ -0,0 +1,20 @@ +=== Test LHAPDF C++ OOP Compatibility Layer === +pid=21, x=0.0001, Q2=100 +Native NeoPDF: 30.8985 +NeoPDF (LHAPDF compat): 30.8985 +Relative difference: 0 +Q2=100 +Native NeoPDF: 0.178123 +NeoPDF (LHAPDF compat): 0.178123 +Relative difference: 0 +LHAPDF C++ OOP Compatibility test passed. +=== Test LHAPDF C Compatibility Layer === +Comparing gluon PDF at x=0.1, Q=10.0 +Native NeoPDF: 1.420505e-01 +NeoPDF (LHAPDF compat): 1.420505e-01 +Relative difference: 0.000000e+00 +Comparing AlphaS at Q=10.0 +Native NeoPDF: 1.781227e-01 +NeoPDF (LHAPDF compat): 1.781227e-01 +Relative difference: 0.000000e+00 +LHAPDF C Compatibility test passed. diff --git a/neopdf_fapi/neopdf.f90 b/neopdf_fapi/neopdf.f90 index 02d2e786..0e2c4853 100644 --- a/neopdf_fapi/neopdf.f90 +++ b/neopdf_fapi/neopdf.f90 @@ -300,6 +300,29 @@ function c_neopdf_grid_compress(collection, metadata, output_path) bind(c, name= character (c_char) :: output_path(*) integer (c_int) :: c_neopdf_grid_compress end function + + subroutine c_setlhaparm(line, len) bind(c, name="setlhaparm_") + use iso_c_binding + character(kind=c_char) :: line(*) + integer(c_int), value :: len + end subroutine + + subroutine c_initpdfsetbyname(name, len) bind(c, name="initpdfsetbyname_") + use iso_c_binding + character(kind=c_char) :: name(*) + integer(c_int), value :: len + end subroutine + + subroutine c_initpdf(member) bind(c, name="initpdf_") + use iso_c_binding + integer(c_int) :: member + end subroutine + + subroutine c_evolvepdf(x, q, f) bind(c, name="evolvepdf_") + use iso_c_binding + real(c_double) :: x, q + real(c_double) :: f(*) + end subroutine end interface contains @@ -586,4 +609,25 @@ function neopdf_grid_compress(collection, metadata, output_path) result(res) res = c_neopdf_grid_compress(collection%ptr, metadata, output_path // c_null_char) end function + subroutine setlhaparm(line) + character(len=*), intent(in) :: line + call c_setlhaparm(line, int(len(trim(line)), c_int)) + end subroutine + + subroutine initpdfsetbyname(name) + character(len=*), intent(in) :: name + call c_initpdfsetbyname(name, int(len(trim(name)), c_int)) + end subroutine + + subroutine initpdf(member) + integer, intent(in) :: member + call c_initpdf(member) + end subroutine + + subroutine evolvepdf(x, q, f) + real(dp), intent(in) :: x, q + real(dp), intent(out) :: f(-6:6) + call c_evolvepdf(x, q, f) + end subroutine + end module