From f08401e7a0f40525340fea9ce88ae220466f25c5 Mon Sep 17 00:00:00 2001 From: Nabin Chapagain Date: Wed, 15 Oct 2025 18:22:48 -0500 Subject: [PATCH] change in getppm, updated cmakelist in xrc folder, new cpp file for pm_sampling --- .../main/+pm/+sampling/@Paradram/getppm.m | 2 +- src/matlab/xrc/CMakeLists.txt | 11 +- src/matlab/xrc/pm_sampling.cpp | 209 ++++++++++++++++++ 3 files changed, 220 insertions(+), 2 deletions(-) create mode 100644 src/matlab/xrc/pm_sampling.cpp diff --git a/src/matlab/main/+pm/+sampling/@Paradram/getppm.m b/src/matlab/main/+pm/+sampling/@Paradram/getppm.m index aa587639..04b64fd4 100644 --- a/src/matlab/main/+pm/+sampling/@Paradram/getppm.m +++ b/src/matlab/main/+pm/+sampling/@Paradram/getppm.m @@ -17,7 +17,7 @@ %> \author %> \AmirShahmoradi, September 1, 2012, 12:00 AM, National Institute for Fusion Studies, The University of Texas at Austin
function ppm = getppm(self) - baselink = self._weblinks.docs.generic.usage.sampling.paradram.quickstart.url + baselink = self.weblinks.docs.generic.usage.sampling.paradram.quickstart.url; ppm = getppm@pm.sampling.Sampler(self) ... + "Use the following object method to read the generated basic output chain file and and unroll the contents as a Markov Chain: " + newline ... + newline ... diff --git a/src/matlab/xrc/CMakeLists.txt b/src/matlab/xrc/CMakeLists.txt index 0d170226..380f6ac1 100644 --- a/src/matlab/xrc/CMakeLists.txt +++ b/src/matlab/xrc/CMakeLists.txt @@ -14,6 +14,7 @@ #################################################################################################################################### #################################################################################################################################### +enable_language(CXX) unset(XRC_FILES) set(XRC_FILES pm_sampling @@ -60,9 +61,17 @@ add_custom_target(matlabmex) unset(xtargets) foreach(xtarget ${XRC_FILES}) set(xlibname "${xtarget}") - set(xrcpath "${CMAKE_CURRENT_SOURCE_DIR}/${xtarget}.c") + set(xrcpath "${CMAKE_CURRENT_SOURCE_DIR}/${xtarget}.cpp") + + # EXPLICITLY SET THE LANGUAGE TO C++ + set_source_files_properties("${xrcpath}" PROPERTIES LANGUAGE CXX) + #set(xlibname "${xtarget}_${build}_${mem}_${parname}_${checking}") matlab_add_mex(NAME "${xtarget}" SHARED SRC "${xrcpath}" OUTPUT_NAME "${xlibname}" LINK_TO "${libname}" R2018a) + + # Ensure the target language is C++ + set_target_properties("${xtarget}" PROPERTIES LINKER_LANGUAGE CXX) + target_compile_definitions("${xtarget}" PUBLIC "OMP_ENABLED=${OMP_ENABLED}") set_target_properties("${xtarget}" PROPERTIES BUILD_WITH_INSTALL_RPATH ON) set_target_properties("${xtarget}" PROPERTIES POSITION_INDEPENDENT_CODE ON) diff --git a/src/matlab/xrc/pm_sampling.cpp b/src/matlab/xrc/pm_sampling.cpp new file mode 100644 index 00000000..d4664d58 --- /dev/null +++ b/src/matlab/xrc/pm_sampling.cpp @@ -0,0 +1,209 @@ +// pm_sampling_mex.cpp — C++ MEX modernization of pm_sampling.c +// +// Matches semantics of the original C MEX file while using the C++ MEX API. +// - Input args (same order/requirements): +// 0) method : char row vector (expects "ParaDRAM") +// 1) fun : function_handle +// 2) ndim : scalar integer +// 3) input : char row vector (passed through to runParaDRAM) +// - No outputs; errors are thrown like the original. +// +// Build: mex -v -R2018a CXXFLAGS="$CXXFLAGS -std=c++17" pm_sampling_mex.cpp + +#include +#include +#include +#include + +#include "mex.hpp" +#include "mexAdapter.hpp" + +using matlab::mex::ArgumentList; +using matlab::data::Array; +using matlab::data::ArrayFactory; +using matlab::data::ArrayType; +using matlab::data::CharArray; +using matlab::data::TypedArray; +using matlab::engine::MATLABEngine; + +/////////////////////////////////////////////////////////////////////////////////////////////////// +// External sampler symbol (name alias kept to mirror the original C file) +/////////////////////////////////////////////////////////////////////////////////////////////////// +#define runParaDRAM runParaDRAMD + +#if OMP_ENABLED +extern "C" int32_t runParaDRAM( double(*getLogFunc)( double logFuncState[] + , int32_t ndim + , int32_t njob + , double *avgTimePerFunCallComp + , double *avgTimePerFunCallComm) + , const int32_t ndim + , const char* input ); +#else +extern "C" int32_t runParaDRAM( double(*getLogFunc)( double state[] + , int32_t ndim ) + , const int32_t ndim + , const char* input ); +#endif + +/////////////////////////////////////////////////////////////////////////////////////////////////// +// MexFunction class +/////////////////////////////////////////////////////////////////////////////////////////////////// +class MexFunction : public matlab::mex::Function { +public: + // Store engine and function handle so static callbacks can use them + static MexFunction* self; + + MexFunction() { self = this; } + ~MexFunction() override = default; + + void operator()(ArgumentList outputs, ArgumentList inputs) override { + engine_ = getEngine(); + ArrayFactory f; + + // === Validate counts === + if (!outputs.empty()) { + err(u"Internal ParaMonte MATLAB library error occurred: Too many output arguments."); + } + if (inputs.size() != 4) { + err(u"Internal ParaMonte MATLAB library error occurred: input variable mismatch."); + } + + // === Parse input #0: method (row char) === + if (inputs[0].getType() != ArrayType::CHAR || inputs[0].getDimensions().size() != 2 || + inputs[0].getDimensions()[0] != 1) { + err(u"Input must be a row vector of characters."); + } + CharArray methodChar = inputs[0]; + const std::u16string method16 = methodChar.toUTF16(); + // We only accept "ParaDRAM" like the original + if (method16.size() < 8 || method16.substr(0,8) != u"ParaDRAM") { + err(u"Internal ParaMonte MATLAB library error occurred: Invalid input sampling method."); + } + + // === Parse input #1: MATLAB function handle === + if (inputs[1].getType() != ArrayType::HANDLE_OBJECT_REF) { + err(u"The second input argument must be a function handle."); + } + funcHandle_ = inputs[1]; // keep a reference-managed handle + + // === Parse input #2: ndim (scalar) === + if (inputs[2].getType() != ArrayType::DOUBLE || inputs[2].getNumberOfElements() != 1) { + err(u"Internal ParaMonte MATLAB library error occurred: Input #2 (ndim) must be a scalar."); + } + const int32_t ndim = static_cast( static_cast>(inputs[2])[0] ); + + // === Parse input #3: input string (row char) === + if (inputs[3].getType() != ArrayType::CHAR || inputs[3].getDimensions().size() != 2 || + inputs[3].getDimensions()[0] != 1) { + err(u"Internal ParaMonte MATLAB library error occurred: Input #3 must be a row char vector."); + } + CharArray inChar = inputs[3]; + inputStr_ = inChar.toAscii(); // ASCII/UTF-8; mirrors mxArrayToString behavior + const char* c_input = inputStr_.c_str(); + + // === Call external sampler === + int32_t stat = 0; +#if OMP_ENABLED + stat = runParaDRAM(&MexFunction::getLogFuncOMP, ndim, c_input); +#else + stat = runParaDRAM(&MexFunction::getLogFuncScalar, ndim, c_input); +#endif + if (stat != 0) { + err(u"Mex:ParaMonte: Runtime Error Occurred."); + } + } + +private: + // Convenience error helper + [[noreturn]] void err(std::u16string msg) { + ArrayFactory f; + engine_->feval(u"error", 0, std::vector{ f.createScalar(msg) }); + throw std::runtime_error("unreachable"); + } + + // === Static callbacks that match runParaDRAM expectations === +#if OMP_ENABLED + // getLogFunc for OMP-enabled multi-job path + static double getLogFuncOMP(double logFuncState[], int32_t ndim, int32_t njob, + double* avgTimePerFunCallComp, double* avgTimePerFunCallComm) + { + // logFuncState layout: for each job j, slot 0 is for logf output, slots 1..ndim are the state + // We need to build an ndim-by-njob matrix of states (without the first column per job). + ArrayFactory f; + auto eng = self->engine_; + + // Construct state matrix (ndim x njob) + matlab::data::TypedArray state = f.createArray({ static_cast(ndim), + static_cast(njob) }); + { + // Fill column-wise like original code + // original loops: for (ijob) for (idim = ndimp1*ijob+1 .. ndimp1*(ijob+1)-1) + const int32_t ndimp1 = ndim + 1; + size_t idx = 0; + for (int32_t j = 0; j < njob; ++j) { + const int32_t base = ndimp1 * j; + for (int32_t d = 1; d < ndimp1; ++d) { + // MATLAB column-major order: state(d, j+1) + // We can set by linear iterator in createArray order (column-major) + state.begin()[idx++] = logFuncState[base + d]; + } + } + } + + // Call feval(funHandle, state) expecting 3 outputs: + // 1) logf (1 x njob double) + // 2) avgTimePerFunCallComp (scalar) + // 3) avgTimePerFunCallComm (scalar) + std::vector in{ self->funcHandle_, state }; + auto out = eng->feval(u"feval", 3, in); + + // Parse outputs + { + // out[0] : vector of length njob (row or column); read linearly + TypedArray logf = out[0]; + // Write back to logFuncState: position 0 in each job’s block + for (int32_t j = 0; j < njob; ++j) { + logFuncState[j * (ndim + 1)] = logf.begin()[static_cast(j)]; + } + } + + *avgTimePerFunCallComp = static_cast>(out[1])[0]; + *avgTimePerFunCallComm = static_cast>(out[2])[0]; + + // The original returned a dummy double (mold = 0). Preserve that. + return 0.0; + } +#else + // getLogFunc for scalar (non-OMP) path + static double getLogFuncScalar(double stateBuf[], int32_t ndim) + { + ArrayFactory f; + auto eng = self->engine_; + + // Build ndim x 1 column vector + TypedArray state = f.createArray({ static_cast(ndim), + static_cast(1) }); + // Copy data (column vector) + { + auto it = state.begin(); + for (int32_t i = 0; i < ndim; ++i) *it++ = stateBuf[i]; + } + + // feval(funHandle, state) → 1 output (scalar logf) + std::vector in{ self->funcHandle_, state }; + auto out = eng->feval(u"feval", 1, in); + + double logf = static_cast>(out[0])[0]; + return logf; + } +#endif + +private: + std::shared_ptr engine_; + Array funcHandle_; // MATLAB function handle (kept alive for callbacks) + std::string inputStr_; // backing storage to keep c_str() valid during sampler call +}; + +// Define the static self pointer +MexFunction* MexFunction::self = nullptr;