diff --git a/MANIFEST.in b/MANIFEST.in index 086d5f71..a4a255ab 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1,5 +1,5 @@ include README.md -recursive-include cpp/RAT * +recursive-include cpp * prune tests prune */__pycache__ global-exclude .git diff --git a/RATapi/utils/plotting.py b/RATapi/utils/plotting.py index 7b8c03e8..e1291cf0 100644 --- a/RATapi/utils/plotting.py +++ b/RATapi/utils/plotting.py @@ -158,8 +158,7 @@ def plot_ref_sld_helper( layers[-1, 1], # Bulk Out data.subRoughs[i], # roughness layer, - len(layer), - 1.0, + 1, ) sld_plot.plot( diff --git a/RATapi/wrappers.py b/RATapi/wrappers.py index c428c254..c0489227 100644 --- a/RATapi/wrappers.py +++ b/RATapi/wrappers.py @@ -66,12 +66,17 @@ def handle(*args): ) return np.array(output, "float").tolist() else: - output, sub_rough = getattr(self.engine, self.function_name)( + matlab_args = [ np.array(args[0], "float"), # params np.array(args[1], "float"), # bulk in np.array(args[2], "float"), # bulk out float(args[3] + 1), # contrast - float(-1 if len(args) < 5 else args[4] + 1), # domain + ] + if len(args) > 4: + matlab_args.append(float(args[4] + 1)) # domain number + + output, sub_rough = getattr(self.engine, self.function_name)( + *matlab_args, nargout=2, ) return np.array(output, "float").tolist(), float(sub_rough) diff --git a/cpp/includes/defines.h b/cpp/includes/defines.h new file mode 100644 index 00000000..a9a59b7f --- /dev/null +++ b/cpp/includes/defines.h @@ -0,0 +1,737 @@ +#ifndef RAT_DEFINES_H +#define RAT_DEFINES_H + +#include +#include +#include + +namespace py = pybind11; + +const std::string docsProgressEventData = R"(The Python binding for the C++ progressEventData struct. +The progress event shows the percentage completion for the calculation. This can be emitted by +the DREAM algorithm only. + +Parameters +---------- +message : str + The title text for the event. +percent : float + The percentage of the calculation completed (as a number between 0 and 1) +)"; + +struct ProgressEventData +{ + std::string message; + double percent; +}; + +const std::string docsPlotEventData = R"(The Python binding for the C++ plotEventData struct. +The plot event data contains intermediate results from the calculation. This can be emitted +by the Simplex and DE algorithms only. + +Parameters +---------- +reflectivity : list + The reflectivity curves for each contrast, with the same range as the data + (``data_range`` in the contrast's ``Data`` object) +shiftedData : list + The data with scalefactors and background corrections applied. +sldProfiles : list + The SLD profiles for each contrast. +resampledLayers : list + If resampling is used, the SLD for each contrast after resampling has been performed. +subRoughs : np.ndarray[np.float] + The substrate roughness values for each contrast. +resample : np.ndarray[np.float] + An array containing whether each contrast was resampled. +dataPresent : np.ndarray[np.float] + Non-zero values indicates if data is present for the contrast. +modelType : str + The model type for the project. +contrastNames : list + The names of all contrasts in the project. +)"; + +struct PlotEventData +{ + py::list reflectivity; + py::list shiftedData; + py::list sldProfiles; + py::list resampledLayers; + py::array_t subRoughs; + py::array_t resample; + py::array_t dataPresent; + std::string modelType; + py::list contrastNames; +}; + + +const std::string docsPredictionIntervals = R"(The Python binding for the C++ predictionIntervals struct. +The Bayesian prediction intervals for 95% and 65% confidence. + +For ``reflectivity`` and ``sld``, each list item is an array +with five rows. The rows represent: + +- 0: the 5th percentile; +- 1: the 35th percentile; +- 2: the mean value of the interval; +- 3: the 65th percentile; +- 4: the 95th percentile. + +Parameters +---------- +reflectivity : list + The prediction interval data for reflectivity of each contrast. +sld : list + The prediction interval data for SLD of each contrast. +sampleChi : np.ndarray[np.float] + The value of sumChi at each point of the Markov chain. +)"; + +struct PredictionIntervals +{ + py::list reflectivity; + py::list sld; + py::array_t sampleChi; +}; + +const std::string docsConfidenceIntervals = R"(The Python binding for the C++ confidenceIntervals struct. +The 65% and 95% confidence intervals for the best fit results. + +Parameters +---------- +percentile95 : np.ndarray[np.float] + The 95% confidence intervals for each fit parameter. +percentile65 : np.ndarray[np.float] + The 65% confidence intervals for each fit parameter. +mean : np.ndarray[np.float] + The mean values for each fit parameter. +)"; + +struct ConfidenceIntervals +{ + py::array_t percentile95; + py::array_t percentile65; + py::array_t mean; +}; + +const std::string docsNestedSamplerOutput = R"(The Python binding for the C++ nestedSamplerOutput struct. +The output information from the Nested Sampler (ns). + +Parameters +---------- +logZ : float + The natural logarithm of the evidence Z for the parameter values. +logZErr : float + The estimated uncertainty in the final value of logZ. +nestSamples : np.ndarray[np.float] + ``NestedSamplerOutput.nestSamples[i, j]`` represents the values + sampled at iteration ``i``, where this value is: + + - the value sampled for parameter ``j``, for ``j`` in ``0:nParams``, + - the minimum log-likelihood for ``j = nParams + 1``. + +postSamples : np.ndarray[np.float] + The posterior values at the points sampled in ``NestedSamplerOutput.nestSamples``. +)"; + +struct NestedSamplerOutput +{ + real_T logZ; + real_T logZErr; + py::array_t nestSamples; + py::array_t postSamples; +}; + +const std::string docsDreamParams = R"(The Python binding for the C++ dreamParams struct. +The parameters used by the inner DREAM algorithm. + +Parameters +---------- +nParams : float + The number of parameters used by the algorithm. +nChains : float + The number of MCMC chains used by the algorithm. +nGenerations : float + The number of DE generations calculated per iteration. +parallel : bool + Whether the algorithm should run chains in parallel. +CPU : float + The number of processor cores used for parallel chains. +jumpProbability : float + A probability range for the size of jumps when performing subspace sampling. +pUnitGamma : float + The probability that the scaling-down factor of jumps will be ignored + and a larger jump will be taken for one iteration. +nCR : float + The number of crossovers performed each iteration. +delta : float + The number of chain mutation pairs proposed each iteration. +steps : float + The number of MCMC steps to perform between conversion checks. +zeta : float + The ergodicity of the algorithm. +outlier : str + What test should be used to detect outliers. +adaptPCR : bool + Whether the crossover probability for differential evolution should be + adapted by the algorithm as it runs. +thinning : float + The thinning rate of each Markov chain (to reduce memory intensity) +epsilon : float + The cutoff threshold for Approximate Bayesian Computation (if used) +ABC : bool + Whether Approximate Bayesian Computation is used. +IO : bool + Whether the algorithm should perform IO writes of the model in parallel. +storeOutput : bool + Whether output model simulations are performed. +R : np.np.ndarray[np.float] + An array where row ``i`` is the list of chains + with which chain ``i`` can mutate. +)"; + +struct DreamParams +{ + real_T nParams; + real_T nChains; + real_T nGenerations; + boolean_T parallel; + real_T CPU; + real_T jumpProbability; + real_T pUnitGamma; + real_T nCR; + real_T delta; + real_T steps; + real_T zeta; + std::string outlier; + boolean_T adaptPCR; + real_T thinning; + real_T epsilon; + boolean_T ABC; + boolean_T IO; + boolean_T storeOutput; + py::array_t R; +}; + +const std::string docsDreamOutput = R"(The Python binding for the C++ DreamOutput struct. +The diagnostic output information from DREAM. + +Parameters +---------- +allChains : np.ndarray[np.float] + An ``nGenerations`` x ``nParams + 2`` x ``nChains`` size array, + where ``chain_k = DreamOutput.allChains[:, :, k]`` + is the data of chain ``k`` in the final iteration; + for generation i of the final iteration, ``chain_k[i, j]`` represents: + + - the sampled value of parameter ``j`` for ``j in 0:nParams``; + - the associated log-prior for those sampled values for ``j = nParams + 1``; + - the associated log-likelihood for those sampled values for ``j = nParams + 2``. + +outlierChains : np.ndarray[np.float] + A two-column array where ``DreamOutput.AR[i, 1]`` is the index of a chain + and ``DreamOutput.AR[i, 0]`` is the length of that chain when it was removed + for being an outlier. +runtime : float + The runtime of the DREAM algorithm in seconds. +iteration : float + The number of iterations performed. +modelOutput : float + Unused. Will always be 0. +AR : np.ndarray[np.float] + A two-column array where ``DreamOutput.AR[i, 0]`` is an iteration number + and ``DreamOutput.AR[i, 1]`` is the average acceptance rate of chain step + proposals for that iteration. +R_stat : np.ndarray[np.float] + An array where ``DreamOutput.R_stat[i, 0]`` is an iteration number and + ``DreamOutput.R_stat[i, j]`` is the convergence statistic for parameter ``j`` + at that iteration (where chains are indexed 1 to ``nParams`` inclusive). +CR : np.ndarray[np.float] + A four-column array where ``DreamOutput.CR[i, 0]`` is an iteration number, + ``and DreamOutput.CR[i, j]`` is the selection probability of the ``j``'th crossover + value for that iteration. +)"; + +struct DreamOutput +{ + py::array_t allChains; + py::array_t outlierChains; + real_T runtime; + real_T iteration; + real_T modelOutput; + py::array_t AR; + py::array_t R_stat; + py::array_t CR; +}; + +const std::string docsBayesResults = R"(The Python binding for the C++ bayesResults struct. +The results of a Bayesian RAT calculation. + +Parameters +---------- +predictionIntervals : RATapi.rat_core.orePredictionIntervals + The prediction intervals. +confidenceIntervals : RATapi.rat_core.ConfidenceIntervals + The 65% and 95% confidence intervals for the best fit results. +dreamParams : RATapi.rat_core.DreamParams + The parameters used by DREAM, if relevant. +dreamOutput : RATapi.rat_core.DreamOutput + The output from DREAM if DREAM was used. +nestedSamplerOutput : RATapi.rat_core.NestedSamplerOutput + The output from nested sampling if ns was used. +chain : np.ndarray + The MCMC chains for each parameter. + The ``i``'th column of the array contains the chain for parameter ``fitNames[i]``. +)"; + +struct BayesResults +{ + PredictionIntervals predictionIntervals; + ConfidenceIntervals confidenceIntervals; + DreamParams dreamParams; + DreamOutput dreamOutput; + NestedSamplerOutput nestedSamplerOutput; + py::array_t chain; +}; + +const std::string docsCalculation = R"(The Python binding for the C++ calculationResult struct. +The goodness of fit from the Abeles calculation. + +Parameters +---------- +chiValues : np.ndarray[np.float] + The chi-squared value for each contrast. +sumChi : float + The sum of the chiValues array. +)"; + +struct Calculation +{ + py::array_t chiValues; + real_T sumChi; +}; + +const std::string docsContrastParams = R"(The Python binding for the C++ contrastParams struct. +The experimental parameters for each contrast. + +Parameters +---------- +scalefactors : np.ndarray[np.float] + The scalefactor values for each contrast. +bulkIn : np.ndarray[np.float] + The bulk in values for each contrast. +bulkOut : np.ndarray[np.float] + The bulk out values for each contrast. +subRoughs : np.ndarray[np.float] + The substrate roughness values for each contrast. +resample : np.ndarray[np.float] + An array containing whether each contrast was resampled. +)"; + +struct ContrastParams +{ + py::array_t scalefactors; + py::array_t bulkIn; + py::array_t bulkOut; + py::array_t subRoughs; + py::array_t resample; +}; + +const std::string docsOutputResult = R"(The C++ result struct of a RAT calculation. + +Parameters +---------- +reflectivity : list + The reflectivity curves for each contrast, with the same range as the data + (``data_range`` in the contrast's ``Data`` object) +simulation : list + The reflectivity curves for each contrast, which can be a wider range to allow extrapolation + (``simulation_range`` in the contrast's ``Data`` object). +shiftedData : list + The data with scalefactors and background corrections applied. +backgrounds : list + The background for each contrast defined over the simulation range. +resolutions : list + The resolution for each contrast defined over the simulation range. +layerSlds : list + The array of layer parameter values for each contrast. +sldProfiles : list + The SLD profiles for each contrast. +resampledLayers : list + If resampling is used, the SLD for each contrast after resampling has been performed. +calculationResults : RATapi.rat_core.Calculation + The chi-squared fit results from the final calculation and fit. +contrastParams : RATapi.rat_core.ContrastParams + The experimental parameters for the contrasts. +fitParams : np.ndarray[np.float] + The best fit value of the parameter with name ``fitNames[i]``. +fitNames : list[str] + The names of the fit parameters, where ``fitNames[i]`` is the name + of the parameter with value given in ``fitParams[i]``. +)"; + +struct OutputResult { + py::list reflectivity; + py::list simulation; + py::list shiftedData; + py::list backgrounds; + py::list resolutions; + py::list layerSlds; + py::list sldProfiles; + py::list resampledLayers; + Calculation calculationResults {}; + ContrastParams contrastParams {}; + py::array_t fitParams; + py::list fitNames; +}; + +const std::string docsLimits = R"(The Python binding for the C++ limit struct which contains +Min and max values for each parameter defined in the project. + +Parameters +---------- +params : np.ndarray[np.float] + Limits for params in the problem definition. +backgroundParams : np.ndarray[np.float] + Limits for backgroundParams in the problem definition. +scalefactors : np.ndarray[np.float] + Limits for scalefactors in the problem definition. +qzshifts : np.ndarray[np.float] + Limits for qzshifts in the problem definition. +bulkIns : np.ndarray[np.float] + Limits for bulkIns in the problem definition. +bulkOuts : np.ndarray[np.float] + Limits for bulkOuts in the problem definition. +resolutionParams : np.ndarray[np.float] + Limits for resolutionParams in the problem definition. +domainRatios : np.ndarray[np.float] + Limits for domainRatios in the problem definition. +)"; + +struct Limits { + py::array_t params; + py::array_t backgroundParams; + py::array_t scalefactors; + py::array_t qzshifts; + py::array_t bulkIns; + py::array_t bulkOuts; + py::array_t resolutionParams; + py::array_t domainRatios; +}; + +const std::string docsNameStore = R"(The Python binding for the C++ names struct which +contains names of all parameters in the project. + +Parameters +---------- +params : list + Names of params in the problem definition. +backgroundParams : list + Names of backgroundParams in the problem definition. +scalefactors : list + Names of scalefactors in the problem definition. +qzshifts : list + Names of qzshifts in the problem definition. +bulkIns : list + Names of bulkIns in the problem definition. +bulkName: list + Names of bulkOuts in the problem definition. +resolutionParams : list + Names of resolutionParams in the problem definition. +domainRatios : list + Names of domainRatios in the problem definition. +)"; + +struct NameStore { + py::list params; + py::list backgroundParams; + py::list scalefactors; + py::list qzshifts; + py::list bulkIns; + py::list bulkOuts; + py::list resolutionParams; + py::list domainRatios; + py::list contrasts; +}; + + +const std::string docsChecks = R"(The Python binding for the C++ checks struct which contains +flags indicating which parameters should be fitted in the project. + +For each attribute, if index ``i`` is non-zero, then parameter ``i`` in that attribute is fitted, e.g. if ``Checks.scalefactors = [0.0, 1.0, 1.0]``, then the second and third scalefactors are fitted and the first is not. + +Parameters +---------- +params : np.ndarray[np.float] + Non-zero values indicates which params is fitted. +backgroundParams : np.ndarray[np.float] + Non-zero values indicates which backgroundParams is fitted. +scalefactors : np.ndarray[np.float] + Non-zero values indicates which scalefactors is fitted. +qzshifts : np.ndarray[np.float] + Non-zero values indicates which qzshifts is fitted. +bulkIns : np.ndarray[np.float] + Non-zero values indicates which bulkIns is fitted. +bulkOuts : np.ndarray[np.float] + Non-zero values indicates which bulkOuts is fitted. +resolutionParams : np.ndarray[np.float] + Non-zero values indicates which resolutionParams is fitted. +domainRatios : np.ndarray[np.float] + Non-zero values indicates which domainRatios is fitted. +)"; + +struct Checks { + py::array_t params; + py::array_t backgroundParams; + py::array_t qzshifts; + py::array_t scalefactors; + py::array_t bulkIns; + py::array_t bulkOuts; + py::array_t resolutionParams; + py::array_t domainRatios; +}; + +const std::string docsProblemDefinition = R"(The Python binding for the C++ problem struct. + +Parameters +---------- +TF : str + The target function for the calculation which can be 'normal' or 'domains'. +resample : np.ndarray[np.float] + If ``resample[i]`` is non-zero, then contrast ``i`` will be resampled. +data : list + Data for each contrast. +dataPresent : np.ndarray[np.float] + If ``dataPresent[i]`` is non-zero, then contrast ``i`` has experimental data. +dataLimits : list + Data limits for each contrast. +simulationLimits : list; + Simulation for each contrast. +oilChiDataPresent : np.ndarray[np.float] + If ``dataPresent[i]`` is non-zero, then contrast ``i`` has oilChi data. This is currently not being used. +numberOfContrasts : int + Number of contrasts. +geometry : str + The geometry to use which can be 'air/substrate' or 'substrate/liquid' +useImaginary : bool + Indicates whether imaginary component is used for the SLD value in layers, i.e. + absorption is set to True for the project. +repeatLayers : list + Information about repeating layers for each contrast. This is currently not being used. +contrastBackgroundParams : list + Indices of backgroundParams used for each contrast +contrastBackgroundTypes : list + Background type for each contrast. +contrastBackgroundActions : list + Background action for each contrast. +contrastQzshifts : np.ndarray[np.float] + Indices of Qzshifts used for each contrast. This is currently not being used. +contrastScalefactors : np.ndarray[np.float] + Indices of scalefactors used for each contrast. +contrastBulkIns : np.ndarray[np.float] + Indices of BulkIns used for each contrast. +contrastBulkOuts : np.ndarray[np.float] + Indices of BulkIns used for each contrast. +contrastResolutionParams : list + Indices of resolutionParams used for each contrast +contrastResolutionTypes : list + Resolution type for each contrast. +backgroundParams : np.ndarray[np.float] + Background parameter values. +qzshifts : np.ndarray[np.float] + Qzshift values. This currently not being used. +scalefactors : np.ndarray[np.float] + Scalefactors values. +bulkIns : np.ndarray[np.float] + BulkIn values. +bulkOuts : np.ndarray[np.float] + BulkOut values. +resolutionParams : np.ndarray[np.float] + Resolution parameter values. +params : np.ndarray[np.float] + Parameter values. +numberOfLayers : int + Number of layers. +contrastLayers : list + Indices of layers added to the model of each contrast. +layersDetails : list + Indices of parameters in each layer. +customFiles : object + Iterable with custom file functions +modelType : str + The layer model type which can be 'standard layers', 'custom layers', or 'custom xy'. +contrastCustomFiles : np.ndarray[np.float] + Indices of CustomFiles used for each domain contrast +contrastDomainRatios : np.ndarray[np.float] + Indices of DomainRatios used for each domain contrast +domainRatios : np.ndarray[np.float] + Domain ratio values +numberOfDomainContrasts : int + Number of domain contrasts. +domainContrastLayers : list + Indices of layers added to the model of each domain contrast. +fitParams : np.ndarray[np.float] + Values of fitted parameters. +otherParams : np.ndarray[np.float] + Values of non-fitted parameters. +fitLimits : np.ndarray[np.float] + Limits of fitted parameters. +otherLimits : np.ndarray[np.float] + Limits of non fitted parameters. +priorNames : list + Parameter names for for all parameters in the problem definition. +priorValues : np.ndarray[np.float] + Prior type, mu, and sigma for all parameters in the problem definition. +names : RATapi.rat_core.NameStore + Names of all parameters. +checks : RATapi.rat_core.Checks + Flags indicating which parameters should be fitted. +)"; + +struct ProblemDefinition { + std::string TF {}; + py::array_t resample; + py::list data; + py::array_t dataPresent; + py::list dataLimits; + py::list simulationLimits; + py::array_t oilChiDataPresent; + real_T numberOfContrasts; + std::string geometry {}; + boolean_T useImaginary {}; + py::list repeatLayers; + py::list contrastBackgroundParams; + py::list contrastBackgroundTypes; + py::list contrastBackgroundActions; + py::array_t contrastQzshifts; + py::array_t contrastScalefactors; + py::array_t contrastBulkIns; + py::array_t contrastBulkOuts; + py::list contrastResolutionParams; + py::list contrastResolutionTypes; + py::array_t backgroundParams; + py::array_t qzshifts; + py::array_t scalefactors; + py::array_t bulkIns; + py::array_t bulkOuts; + py::array_t resolutionParams; + py::array_t params; + real_T numberOfLayers {}; + py::list contrastLayers; + py::list layersDetails; + py::object customFiles; + std::string modelType {}; + py::array_t contrastCustomFiles; + py::array_t contrastDomainRatios; + py::array_t domainRatios; + real_T numberOfDomainContrasts {}; + py::list domainContrastLayers; + py::array_t fitParams; + py::array_t otherParams; + py::array_t fitLimits; + py::array_t otherLimits; + py::list priorNames; + py::array_t priorValues; + NameStore names; + Checks checks {}; +}; + + +const std::string docsControl = R"(The Python binding for the C++ controls struct. + +Parameters +---------- +parallel : str + How the calculation should be parallelised (This uses the Parallel Computing Toolbox). Can be 'single', 'contrasts' or 'points'. +procedure : str + Which procedure RAT should execute. Can be 'calculate', 'simplex', 'de', 'ns', or 'dream'. +calcSldDuringFit : bool + Whether SLD will be calculated during fit (for live plotting etc.) +resampleMinAngle : float + The upper threshold on the angle between three sampled points for resampling, in units of radians over pi. +resampleNPoints : int + The number of initial points to use for resampling. +display : str + How much RAT should print to the terminal. Can be 'off', 'iter', 'notify', or 'final'. +IPCFilePath : str + The path of the inter process communication file. +updateFreq : int + [SIMPLEX, DE] Number of iterations between printing progress updates to the terminal. +updatePlotFreq : int + [SIMPLEX, DE] Number of iterations between updates to live plots. +xTolerance : float + [SIMPLEX] The termination tolerance for step size. +funcTolerance : float + [SIMPLEX] The termination tolerance for change in chi-squared. +maxFuncEvals : int + [SIMPLEX] The maximum number of function evaluations before the algorithm terminates. +maxIterations : int + [SIMPLEX] The maximum number of iterations before the algorithm terminates. +populationSize : int + [DE] The number of candidate solutions that exist at any time. +fWeight : float + [DE] The step size for how different mutations are to their parents. +crossoverProbability : float + [DE] The probability of exchange of parameters between individuals at any iteration. +strategy : int + [DE] The algorithm used to generate new candidates. +targetValue : float + [DE] The value of chi-squared at which the algorithm will terminate. +numGenerations : int + [DE] The maximum number of iterations before the algorithm terminates. +nLive : int + [NS] The number of points to sample. +nMCMC : int + [NS] If non-zero, an MCMC process with ``nMCMC`` chains will be used instead of MultiNest. +propScale : float + [NS] A scaling factor for the ellipsoid generated by MultiNest. +nsTolerance : float + [NS] The tolerance threshold for when the algorithm should terminate. +nSamples : int + [DREAM] The number of samples in the initial population for each chain. +nChains : int + [DREAM] The number of Markov chains to use in the algorithm. +jumpProbability : float + [DREAM] The probability range for the size of jumps in sampling. Larger values mean more variable jumps. +pUnitGamma : float + [DREAM] The probability that the scaling-down factor of jumps will be ignored and a larger jump will be taken. +boundHandling : str + [DREAM] How steps past the space boundaries should be handled. Can be 'off', 'reflect', 'bound', or 'fold'. +adaptPCR : bool + [DREAM] Whether the crossover probability for differential evolution should be adapted during the run. +)"; + +struct Control { + std::string parallel {}; + std::string procedure {}; + std::string display {}; + real_T xTolerance {}; + real_T funcTolerance {}; + real_T maxFuncEvals {}; + real_T maxIterations {}; + real_T populationSize {}; + real_T fWeight {}; + real_T crossoverProbability {}; + real_T targetValue {}; + real_T numGenerations {}; + real_T strategy {}; + real_T nLive {}; + real_T nMCMC {}; + real_T propScale {}; + real_T nsTolerance {}; + boolean_T calcSldDuringFit {}; + real_T resampleMinAngle {}; + real_T resampleNPoints {}; + real_T updateFreq {}; + real_T updatePlotFreq {}; + real_T nSamples {}; + real_T nChains {}; + real_T jumpProbability {}; + real_T pUnitGamma {}; + std::string boundHandling {}; + boolean_T adaptPCR; + std::string IPCFilePath {}; +}; + +#endif // RAT_DEFINES_H \ No newline at end of file diff --git a/cpp/includes/functions.h b/cpp/includes/functions.h new file mode 100644 index 00000000..3ced72fe --- /dev/null +++ b/cpp/includes/functions.h @@ -0,0 +1,500 @@ +#ifndef RAT_FUNCTIONS_H +#define RAT_FUNCTIONS_H + +#include +#include +#include +#include +#include "../RAT/classHandle.hpp" + +namespace py = pybind11; + +template +auto customCaller(std::string identifier, Function f, Args&& ... args) -> decltype((*f)(std::forward(args)...)) +{ + try + { + return (*f)(std::forward(args)...); + } + catch(const std::runtime_error& re) + { + std::string errorMsg; + size_t start_pos = std::string(re.what()).find("$id"); + if(start_pos == std::string::npos) + { + errorMsg = std::string("Error occurred when setting ") + identifier + ". " + re.what(); + } + else + { + errorMsg = re.what(); + errorMsg.replace(start_pos, 3, identifier); + } + + throw std::runtime_error(errorMsg); + } +} + +class Library: public CallbackInterface +{ + public: + + py::function function; + + Library(const py::function function){ + this->function = function; + }; + + + void setOutput(py::tuple& result, std::vector& output, double *outputSize) + { + int nRows = 0, idx = 0; + for (py::handle rowHandle : result[0]) + { + py::list rows = py::cast(rowHandle); + for (py::handle value : rows) + { + output.push_back(py::cast(value)); + idx++; + } + nRows++; + } + + outputSize[0] = nRows; + outputSize[1] = (nRows == 0) ? 0 : idx / nRows; + } + + // Backgrounds + void invoke(std::vector& xdata, std::vector& params, std::vector& output) + { + auto f = py::cast>(this->function); + auto result = f(py::cast(xdata), py::cast(params)); + for (py::handle rowHandle : result) + { + if (py::isinstance(rowHandle)) + output.push_back(py::cast(py::cast(rowHandle)[0])); + else + output.push_back(py::cast(rowHandle)); + + } + }; + + // Domain overload + void invoke(std::vector& params, std::vector& bulkIn, std::vector& bulkOut, + int contrast, int domainNumber, std::vector& output, double *outputSize, double *roughness) + { + auto f = py::cast>(this->function); + auto result = f(py::cast(params), py::cast(bulkIn), py::cast(bulkOut), contrast, domainNumber); + *roughness = py::cast(result[1]); + setOutput(result, output, outputSize); + }; + + // Non-Domain overload + void invoke(std::vector& params, std::vector& bulkIn, std::vector& bulkOut, + int contrast, std::vector& output, double *outputSize, double *roughness) + { + auto f = py::cast>(this->function); + auto result = f(py::cast(params), py::cast(bulkIn), py::cast(bulkOut), contrast); + *roughness = py::cast(result[1]); + setOutput(result, output, outputSize); + }; +}; + + +void stringToRatBoundedArray(std::string value, char_T result_data[], int32_T result_size[2]) +{ + result_size[0] = 1; + result_size[1] = value.length(); + + for (int32_T idx1{0}; idx1 < value.length(); idx1++) { + result_data[idx1] = value[idx1]; + } +} + +void stringToRatCharArray(std::string value, coder::array& result) +{ + result.set_size(1, value.length()); + + for (int32_T idx{0}; idx < value.length(); idx++) { + result[idx] = value[idx]; + } +} + +void stringFromRatBoundedArray(const char_T array_data[], const int32_T array_size[2], std::string& result) +{ + result.resize(array_size[1]); + memcpy(&result[0], array_data, array_size[1]); +} + + +coder::array pyArrayToRatRowArray1d(py::array_t value) +{ + coder::array result; + + py::buffer_info buffer_info = value.request(); + + if (buffer_info.size == 0) + return result; + + if (buffer_info.ndim != 1) + throw std::runtime_error("Expects a 1D numeric array"); + + result.set_size(1, buffer_info.shape[0]); + for (int32_T idx0{0}; idx0 < buffer_info.shape[0]; idx0++) { + result[idx0] = value.at(idx0); + } + + return result; +} + +coder::bounded_array pyArrayToRatBoundedArray(py::array_t value) +{ + coder::bounded_array result {}; + + py::buffer_info buffer_info = value.request(); + + if (buffer_info.size == 0) + return result; + + if (buffer_info.ndim != 1) + throw std::runtime_error("Expects a 1D numeric array"); + + result.size[0] = 1; + result.size[1] = buffer_info.shape[0]; + for (int32_T idx0{0}; idx0 < buffer_info.shape[0]; idx0++) { + result.data[idx0] = value.at(idx0); + } + + return result; +} + +coder::bounded_array pyArrayToRatBoundedArray3(py::array_t value) +{ + coder::bounded_array result {}; + + py::buffer_info buffer_info = value.request(); + + if (buffer_info.size == 0) + return result; + + if (buffer_info.ndim != 1) + throw std::runtime_error("Expects a 1D numeric array"); + + result.size[0] = 1; + result.size[1] = buffer_info.shape[0]; + for (int32_T idx0{0}; idx0 < buffer_info.shape[0]; idx0++) { + result.data[idx0] = value.at(idx0); + } + + return result; +} + +coder::array pyArrayToRatArray2d(py::array_t value) +{ + coder::array result; + + py::buffer_info buffer_info = value.request(); + + if (buffer_info.size == 0) + return result; + + if (buffer_info.ndim != 2) + throw std::runtime_error("Expects a 2D numeric array"); + + result.set_size(buffer_info.shape[0], buffer_info.shape[1]); + + int32_T idx {0}; + for (int32_T idx0{0}; idx0 < buffer_info.shape[0]; idx0++) { + for (int32_T idx1{0}; idx1 < buffer_info.shape[1]; idx1++) { + idx = idx0 + result.size(0) * idx1; + result[idx] = value.at(idx0, idx1); + } + } + + return result; +} + +coder::array pyListToRatCellWrap1(py::list values) +{ + coder::array result; + result.set_size(1, values.size()); + int32_T idx {0}; + for (py::handle array: values) + { + py::array_t casted_array = py::cast(array); + result[idx].f1 = customCaller("$id[" + std::to_string(idx) +"]", pyArrayToRatArray2d, casted_array); + idx++; + } + + return result; +} + +coder::array pyListToRatCellWrap2(py::list values) +{ + coder::array result; + result.set_size(1, values.size()); + int32_T idx {0}; + for (py::handle array: values) + { + py::array_t casted_array = py::cast(array); + if (casted_array.size() != 2) + throw std::runtime_error("Expects a 2D list where each row contains exactly 2 numbers"); + result[idx].f1[0] = casted_array.at(0); + result[idx].f1[1] = casted_array.at(1); + idx++; + } + + return result; +} + + +coder::array pyListToRatCellWrap3(py::list values) +{ + coder::array result; + result.set_size(1, values.size()); + int32_T idx {0}; + for (py::handle array: values) + { + py::array_t casted_array = py::cast(array); + result[idx].f1 = customCaller("$id[" + std::to_string(idx) +"]", pyArrayToRatBoundedArray3, casted_array); + idx++; + } + + return result; +} + +coder::array pyListToRatCellWrap4(py::list values) +{ + coder::array result; + result.set_size(1, values.size()); + int32_T idx {0}; + for (py::handle array: values) + { + py::array_t casted_array = py::cast(array); + result[idx].f1 = customCaller("$id[" + std::to_string(idx) +"]", pyArrayToRatBoundedArray3, casted_array); + idx++; + } + + return result; +} + +coder::array pyListToRatCellWrap5(py::list values) +{ + coder::array result; + result.set_size(1, values.size()); + int32_T idx {0}; + for (py::handle array: values) + { + py::array_t casted_array = py::cast(array); + result[idx].f1 = customCaller("$id[" + std::to_string(idx) +"]", pyArrayToRatRowArray1d, casted_array); + idx++; + } + + return result; +} + +coder::array pyListToRatCellWrap6(py::list values) +{ + coder::array result; + result.set_size(1, values.size()); + int32_T idx {0}; + for (py::handle array: values) + { + py::array_t casted_array = py::cast(array); + result[idx].f1 = customCaller("$id[" + std::to_string(idx) +"]", pyArrayToRatBoundedArray, casted_array); + idx++; + } + + return result; +} + +coder::array pyListToRatCellWrap01d(py::list values) +{ + coder::array result; + result.set_size(values.size()); + int32_T idx {0}; + for (py::handle array: values) + { + if (py::isinstance(array)) { + std::string name = py::cast(array); + stringToRatBoundedArray(name, result[idx].f1.data, result[idx].f1.size); + idx++; + } + else + throw std::runtime_error("Expects a 1D list of strings"); + } + + return result; +} + +coder::array pyListToRatCellWrap02d(py::list values) +{ + coder::array result; + result.set_size(1, values.size()); + int32_T idx {0}; + for (py::handle array: values) + { + if (py::isinstance(array)) { + std::string name = py::cast(array); + stringToRatBoundedArray(name, result[idx].f1.data, result[idx].f1.size); + idx++; + } + else + throw std::runtime_error("Expects a 1D list of strings"); + } + + return result; +} + +coder::array py_function_array_to_rat_cell_wrap_0(py::object values) +{ + auto handles = py::cast(values); + coder::array result; + result.set_size(1, handles.size()); + int32_T idx {0}; + for (py::handle array: handles) + { + auto func = py::cast(array); + std::string func_ptr = convertPtr2String(new Library(func)); + stringToRatBoundedArray(func_ptr, result[idx].f1.data, result[idx].f1.size); + idx++; + } + + return result; +} +template +py::array_t pyArrayFromRatArray1d(T array, bool isCol=true) +{ + auto size = isCol ? array.size(1) : array.size(0); + auto result_array = py::array_t(size); + std::memcpy(result_array.request().ptr, array.data(), result_array.nbytes()); + + return result_array; +} + +py::array_t pyArrayFromRatArray2d(coder::array array) +{ + auto result_array = py::array_t({array.size(0), array.size(1)}); + std::memcpy(result_array.request().ptr, array.data(), result_array.nbytes()); + + return result_array; +} + +py::list pyListFromRatCellWrap01d(coder::array values) +{ + py::list result; + for (int32_T idx0{0}; idx0 < values.size(0); idx0++) { + std::string tmp; + stringFromRatBoundedArray(values[idx0].f1.data, values[idx0].f1.size, tmp); + result.append(tmp); + } + + return result; +} + +py::list pyListFromRatCellWrap02d(coder::array values) +{ + py::list result; + for (int32_T idx0{0}; idx0 < values.size(1); idx0++) { + std::string tmp; + stringFromRatBoundedArray(values[idx0].f1.data, values[idx0].f1.size, tmp); + result.append(tmp); + } + + return result; +} + +py::list pyListFromRatCellWrap2(coder::array values) +{ + py::list result; + + for (int32_T idx0{0}; idx0 < values.size(1); idx0++) { + py::list inner = py::make_tuple(values[idx0].f1[0], values[idx0].f1[1]); + result.append(inner); + } + + return result; +} + +template +py::list pyList1DFromRatCellWrap2D(const T& values) +{ + py::list result; + + for (int32_T idx0{0}; idx0 < values.size(1); idx0++) { + result.append(pyArrayFromRatArray2d(values[idx0].f1)); + } + + return result; +} + +template +py::list pyList1DFromRatCellWrap1D(const T& values) +{ + py::list result; + + for (int32_T idx0{0}; idx0 < values.size(0); idx0++) { + result.append(pyArrayFromRatArray2d(values[idx0].f1)); + } + + return result; +} + +template +py::list pyList2dFromRatCellWrap(const T& values) +{ + py::list result; + int32_T idx {0}; + for (int32_T idx0{0}; idx0 < values.size(0); idx0++) { + py::list inner; + for (int32_T idx1{0}; idx1 < values.size(1); idx1++) { + idx = idx0 + values.size(0) * idx1; + inner.append(pyArrayFromRatArray2d(values[idx].f1)); + } + result.append(inner); + } + + return result; +} + +template +py::list pyListFromBoundedCellWrap(const T& values) +{ + py::list result; + + for (int32_T idx0{0}; idx0 < values.size(1); idx0++) { + auto array = py::array_t({values[idx0].f1.size[0]}); + std::memcpy(array.request().ptr, values[idx0].f1.data, array.nbytes()); + + result.append(array); + } + + return result; +} + +template +py::array_t pyArray1dFromBoundedArray(const T& array) +{ + auto result_array = py::array_t({array.size[0]}); + std::memcpy(result_array.request().ptr, array.data, result_array.nbytes()); + + return result_array; +} + +template +py::array_t pyArray2dFromBoundedArray(const T& array) +{ + auto result_array = py::array_t({array.size[0], array.size[1]}); + std::memcpy(result_array.request().ptr, array.data, result_array.nbytes()); + + return result_array; +} + +py::array_t pyArrayFromRatArray3d(coder::array array) +{ + auto result_array = py::array_t({array.size(0), array.size(1), array.size(2)}); + std::memcpy(result_array.request().ptr, array.data(), result_array.nbytes()); + + return result_array; +} + +#endif // RAT_FUNCTIONS_H \ No newline at end of file diff --git a/cpp/rat.cpp b/cpp/rat.cpp index f159927f..84c30243 100644 --- a/cpp/rat.cpp +++ b/cpp/rat.cpp @@ -17,103 +17,15 @@ setup_pybind11(cfg) #include "RAT/RATMain_terminate.h" #include "RAT/RATMain_types.h" #include "RAT/makeSLDProfileXY.h" -#include "RAT/classHandle.hpp" #include "RAT/dylib.hpp" #include "RAT/events/eventManager.h" +#include "includes/defines.h" +#include "includes/functions.h" namespace py = pybind11; const int DEFAULT_DOMAIN = -1; - -template -auto customCaller(std::string identifier, Function f, Args&& ... args) -> decltype((*f)(std::forward(args)...)) -{ - try - { - return (*f)(std::forward(args)...); - } - catch(const std::runtime_error& re) - { - std::string errorMsg; - size_t start_pos = std::string(re.what()).find("$id"); - if(start_pos == std::string::npos) - { - errorMsg = std::string("Error occurred when setting ") + identifier + ". " + re.what(); - } - else - { - errorMsg = re.what(); - errorMsg.replace(start_pos, 3, identifier); - } - - throw std::runtime_error(errorMsg); - } -} - -class Library: public CallbackInterface -{ - public: - - py::function function; - - Library(const py::function function){ - this->function = function; - }; - - - void setOutput(py::tuple& result, std::vector& output, double *outputSize) - { - int nRows = 0, idx = 0; - for (py::handle rowHandle : result[0]) - { - py::list rows = py::cast(rowHandle); - for (py::handle value : rows) - { - output.push_back(py::cast(value)); - idx++; - } - nRows++; - } - - outputSize[0] = nRows; - outputSize[1] = (nRows == 0) ? 0 : idx / nRows; - } - - // Backgrounds - void invoke(std::vector& xdata, std::vector& params, std::vector& output) - { - auto f = py::cast>(this->function); - auto result = f(py::cast(xdata), py::cast(params)); - for (py::handle rowHandle : result) - { - if (py::isinstance(rowHandle)) - output.push_back(py::cast(py::cast(rowHandle)[0])); - else - output.push_back(py::cast(rowHandle)); - - } - }; - - // Domain overload - void invoke(std::vector& params, std::vector& bulkIn, std::vector& bulkOut, - int contrast, int domainNumber, std::vector& output, double *outputSize, double *roughness) - { - auto f = py::cast>(this->function); - auto result = f(py::cast(params), py::cast(bulkIn), py::cast(bulkOut), contrast, domainNumber); - *roughness = py::cast(result[1]); - setOutput(result, output, outputSize); - }; - - // Non-Domain overload - void invoke(std::vector& params, std::vector& bulkIn, std::vector& bulkOut, - int contrast, std::vector& output, double *outputSize, double *roughness) - { - auto f = py::cast>(this->function); - auto result = f(py::cast(params), py::cast(bulkIn), py::cast(bulkOut), contrast); - *roughness = py::cast(result[1]); - setOutput(result, output, outputSize); - }; -}; +const int DEFAULT_NREPEATS = 1; class DylibEngine { @@ -185,25 +97,6 @@ class DylibEngine }; }; -struct ProgressEventData -{ - std::string message; - double percent; -}; - -struct PlotEventData -{ - py::list reflectivity; - py::list shiftedData; - py::list sldProfiles; - py::list resampledLayers; - py::array_t subRoughs; - py::array_t resample; - py::array_t dataPresent; - std::string modelType; - py::list contrastNames; -}; - class EventBridge { public: @@ -328,480 +221,6 @@ class EventBridge }; }; -struct PredictionIntervals -{ - py::list reflectivity; - py::list sld; - py::array_t sampleChi; -}; - -struct ConfidenceIntervals -{ - py::array_t percentile95; - py::array_t percentile65; - py::array_t mean; -}; - -struct NestedSamplerOutput -{ - real_T logZ; - real_T logZErr; - py::array_t nestSamples; - py::array_t postSamples; -}; - -struct DreamParams -{ - real_T nParams; - real_T nChains; - real_T nGenerations; - boolean_T parallel; - real_T CPU; - real_T jumpProbability; - real_T pUnitGamma; - real_T nCR; - real_T delta; - real_T steps; - real_T zeta; - std::string outlier; - boolean_T adaptPCR; - real_T thinning; - real_T epsilon; - boolean_T ABC; - boolean_T IO; - boolean_T storeOutput; - py::array_t R; -}; - -struct DreamOutput -{ - py::array_t allChains; - py::array_t outlierChains; - real_T runtime; - real_T iteration; - real_T modelOutput; - py::array_t AR; - py::array_t R_stat; - py::array_t CR; -}; - -struct BayesResults -{ - PredictionIntervals predictionIntervals; - ConfidenceIntervals confidenceIntervals; - DreamParams dreamParams; - DreamOutput dreamOutput; - NestedSamplerOutput nestedSamplerOutput; - py::array_t chain; -}; - -struct Checks { - py::array_t params; - py::array_t backgroundParams; - py::array_t qzshifts; - py::array_t scalefactors; - py::array_t bulkIns; - py::array_t bulkOuts; - py::array_t resolutionParams; - py::array_t domainRatios; -}; - -struct Calculation -{ - py::array_t chiValues; - real_T sumChi; -}; - -struct ContrastParams -{ - py::array_t scalefactors; - py::array_t bulkIn; - py::array_t bulkOut; - py::array_t subRoughs; - py::array_t resample; -}; - -struct OutputResult { - py::list reflectivity; - py::list simulation; - py::list shiftedData; - py::list backgrounds; - py::list resolutions; - py::list layerSlds; - py::list sldProfiles; - py::list resampledLayers; - Calculation calculationResults {}; - ContrastParams contrastParams {}; - py::array_t fitParams; - py::list fitNames; -}; - -struct Limits { - py::array_t params; - py::array_t backgroundParams; - py::array_t scalefactors; - py::array_t qzshifts; - py::array_t bulkIns; - py::array_t bulkOuts; - py::array_t resolutionParams; - py::array_t domainRatios; -}; - -struct NameStore { - py::list params; - py::list backgroundParams; - py::list scalefactors; - py::list qzshifts; - py::list bulkIns; - py::list bulkOuts; - py::list resolutionParams; - py::list domainRatios; - py::list contrasts; -}; - -struct ProblemDefinition { - std::string TF {}; - py::array_t resample; - py::list data; - py::array_t dataPresent; - py::list dataLimits; - py::list simulationLimits; - py::array_t oilChiDataPresent; - real_T numberOfContrasts; - std::string geometry {}; - boolean_T useImaginary {}; - py::list repeatLayers; - py::list contrastBackgroundParams; - py::list contrastBackgroundTypes; - py::list contrastBackgroundActions; - py::array_t contrastQzshifts; - py::array_t contrastScalefactors; - py::array_t contrastBulkIns; - py::array_t contrastBulkOuts; - py::list contrastResolutionParams; - py::list contrastResolutionTypes; - py::array_t backgroundParams; - py::array_t qzshifts; - py::array_t scalefactors; - py::array_t bulkIns; - py::array_t bulkOuts; - py::array_t resolutionParams; - py::array_t params; - real_T numberOfLayers {}; - py::list contrastLayers; - py::list layersDetails; - py::object customFiles; - std::string modelType {}; - py::array_t contrastCustomFiles; - py::array_t contrastDomainRatios; - py::array_t domainRatios; - real_T numberOfDomainContrasts {}; - py::list domainContrastLayers; - py::array_t fitParams; - py::array_t otherParams; - py::array_t fitLimits; - py::array_t otherLimits; - py::list priorNames; - py::array_t priorValues; - NameStore names; - Checks checks {}; -}; - -struct Control { - std::string parallel {}; - std::string procedure {}; - std::string display {}; - real_T xTolerance {}; - real_T funcTolerance {}; - real_T maxFuncEvals {}; - real_T maxIterations {}; - real_T populationSize {}; - real_T fWeight {}; - real_T crossoverProbability {}; - real_T targetValue {}; - real_T numGenerations {}; - real_T strategy {}; - real_T nLive {}; - real_T nMCMC {}; - real_T propScale {}; - real_T nsTolerance {}; - boolean_T calcSldDuringFit {}; - real_T resampleMinAngle {}; - real_T resampleNPoints {}; - real_T updateFreq {}; - real_T updatePlotFreq {}; - real_T nSamples {}; - real_T nChains {}; - real_T jumpProbability {}; - real_T pUnitGamma {}; - std::string boundHandling {}; - boolean_T adaptPCR; - std::string IPCFilePath {}; -}; - - -void stringToRatBoundedArray(std::string value, char_T result_data[], int32_T result_size[2]) -{ - result_size[0] = 1; - result_size[1] = value.length(); - - for (int32_T idx1{0}; idx1 < value.length(); idx1++) { - result_data[idx1] = value[idx1]; - } -} - -void stringToRatCharArray(std::string value, coder::array& result) -{ - result.set_size(1, value.length()); - - for (int32_T idx{0}; idx < value.length(); idx++) { - result[idx] = value[idx]; - } -} - -void stringFromRatBoundedArray(const char_T array_data[], const int32_T array_size[2], std::string& result) -{ - result.resize(array_size[1]); - memcpy(&result[0], array_data, array_size[1]); -} - - -coder::array pyArrayToRatRowArray1d(py::array_t value) -{ - coder::array result; - - py::buffer_info buffer_info = value.request(); - - if (buffer_info.size == 0) - return result; - - if (buffer_info.ndim != 1) - throw std::runtime_error("Expects a 1D numeric array"); - - result.set_size(1, buffer_info.shape[0]); - for (int32_T idx0{0}; idx0 < buffer_info.shape[0]; idx0++) { - result[idx0] = value.at(idx0); - } - - return result; -} - -coder::bounded_array pyArrayToRatBoundedArray(py::array_t value) -{ - coder::bounded_array result {}; - - py::buffer_info buffer_info = value.request(); - - if (buffer_info.size == 0) - return result; - - if (buffer_info.ndim != 1) - throw std::runtime_error("Expects a 1D numeric array"); - - result.size[0] = 1; - result.size[1] = buffer_info.shape[0]; - for (int32_T idx0{0}; idx0 < buffer_info.shape[0]; idx0++) { - result.data[idx0] = value.at(idx0); - } - - return result; -} - -coder::bounded_array pyArrayToRatBoundedArray3(py::array_t value) -{ - coder::bounded_array result {}; - - py::buffer_info buffer_info = value.request(); - - if (buffer_info.size == 0) - return result; - - if (buffer_info.ndim != 1) - throw std::runtime_error("Expects a 1D numeric array"); - - result.size[0] = 1; - result.size[1] = buffer_info.shape[0]; - for (int32_T idx0{0}; idx0 < buffer_info.shape[0]; idx0++) { - result.data[idx0] = value.at(idx0); - } - - return result; -} - -coder::array pyArrayToRatArray2d(py::array_t value) -{ - coder::array result; - - py::buffer_info buffer_info = value.request(); - - if (buffer_info.size == 0) - return result; - - if (buffer_info.ndim != 2) - throw std::runtime_error("Expects a 2D numeric array"); - - result.set_size(buffer_info.shape[0], buffer_info.shape[1]); - - int32_T idx {0}; - for (int32_T idx0{0}; idx0 < buffer_info.shape[0]; idx0++) { - for (int32_T idx1{0}; idx1 < buffer_info.shape[1]; idx1++) { - idx = idx0 + result.size(0) * idx1; - result[idx] = value.at(idx0, idx1); - } - } - - return result; -} - -coder::array pyListToRatCellWrap1(py::list values) -{ - coder::array result; - result.set_size(1, values.size()); - int32_T idx {0}; - for (py::handle array: values) - { - py::array_t casted_array = py::cast(array); - result[idx].f1 = customCaller("$id[" + std::to_string(idx) +"]", pyArrayToRatArray2d, casted_array); - idx++; - } - - return result; -} - -coder::array pyListToRatCellWrap2(py::list values) -{ - coder::array result; - result.set_size(1, values.size()); - int32_T idx {0}; - for (py::handle array: values) - { - py::array_t casted_array = py::cast(array); - if (casted_array.size() != 2) - throw std::runtime_error("Expects a 2D list where each row contains exactly 2 numbers"); - result[idx].f1[0] = casted_array.at(0); - result[idx].f1[1] = casted_array.at(1); - idx++; - } - - return result; -} - - -coder::array pyListToRatCellWrap3(py::list values) -{ - coder::array result; - result.set_size(1, values.size()); - int32_T idx {0}; - for (py::handle array: values) - { - py::array_t casted_array = py::cast(array); - result[idx].f1 = customCaller("$id[" + std::to_string(idx) +"]", pyArrayToRatBoundedArray3, casted_array); - idx++; - } - - return result; -} - -coder::array pyListToRatCellWrap4(py::list values) -{ - coder::array result; - result.set_size(1, values.size()); - int32_T idx {0}; - for (py::handle array: values) - { - py::array_t casted_array = py::cast(array); - result[idx].f1 = customCaller("$id[" + std::to_string(idx) +"]", pyArrayToRatBoundedArray3, casted_array); - idx++; - } - - return result; -} - -coder::array pyListToRatCellWrap5(py::list values) -{ - coder::array result; - result.set_size(1, values.size()); - int32_T idx {0}; - for (py::handle array: values) - { - py::array_t casted_array = py::cast(array); - result[idx].f1 = customCaller("$id[" + std::to_string(idx) +"]", pyArrayToRatRowArray1d, casted_array); - idx++; - } - - return result; -} - -coder::array pyListToRatCellWrap6(py::list values) -{ - coder::array result; - result.set_size(1, values.size()); - int32_T idx {0}; - for (py::handle array: values) - { - py::array_t casted_array = py::cast(array); - result[idx].f1 = customCaller("$id[" + std::to_string(idx) +"]", pyArrayToRatBoundedArray, casted_array); - idx++; - } - - return result; -} - -coder::array pyListToRatCellWrap01d(py::list values) -{ - coder::array result; - result.set_size(values.size()); - int32_T idx {0}; - for (py::handle array: values) - { - if (py::isinstance(array)) { - std::string name = py::cast(array); - stringToRatBoundedArray(name, result[idx].f1.data, result[idx].f1.size); - idx++; - } - else - throw std::runtime_error("Expects a 1D list of strings"); - } - - return result; -} - -coder::array pyListToRatCellWrap02d(py::list values) -{ - coder::array result; - result.set_size(1, values.size()); - int32_T idx {0}; - for (py::handle array: values) - { - if (py::isinstance(array)) { - std::string name = py::cast(array); - stringToRatBoundedArray(name, result[idx].f1.data, result[idx].f1.size); - idx++; - } - else - throw std::runtime_error("Expects a 1D list of strings"); - } - - return result; -} - -coder::array py_function_array_to_rat_cell_wrap_0(py::object values) -{ - auto handles = py::cast(values); - coder::array result; - result.set_size(1, handles.size()); - int32_T idx {0}; - for (py::handle array: handles) - { - auto func = py::cast(array); - std::string func_ptr = convertPtr2String(new Library(func)); - stringToRatBoundedArray(func_ptr, result[idx].f1.data, result[idx].f1.size); - idx++; - } - - return result; -} - RAT::struct1_T createStruct1(const NameStore& names) { RAT::struct1_T names_struct; @@ -941,143 +360,6 @@ RAT::struct4_T createStruct4(const Control& control) return control_struct; } - -template -py::array_t pyArrayFromRatArray1d(T array, bool isCol=true) -{ - auto size = isCol ? array.size(1) : array.size(0); - auto result_array = py::array_t(size); - std::memcpy(result_array.request().ptr, array.data(), result_array.nbytes()); - - return result_array; -} - -py::array_t pyArrayFromRatArray2d(coder::array array) -{ - auto result_array = py::array_t({array.size(0), array.size(1)}); - std::memcpy(result_array.request().ptr, array.data(), result_array.nbytes()); - - return result_array; -} - -py::list pyListFromRatCellWrap01d(coder::array values) -{ - py::list result; - for (int32_T idx0{0}; idx0 < values.size(0); idx0++) { - std::string tmp; - stringFromRatBoundedArray(values[idx0].f1.data, values[idx0].f1.size, tmp); - result.append(tmp); - } - - return result; -} - -py::list pyListFromRatCellWrap02d(coder::array values) -{ - py::list result; - for (int32_T idx0{0}; idx0 < values.size(1); idx0++) { - std::string tmp; - stringFromRatBoundedArray(values[idx0].f1.data, values[idx0].f1.size, tmp); - result.append(tmp); - } - - return result; -} - -py::list pyListFromRatCellWrap2(coder::array values) -{ - py::list result; - - for (int32_T idx0{0}; idx0 < values.size(1); idx0++) { - py::list inner = py::make_tuple(values[idx0].f1[0], values[idx0].f1[1]); - result.append(inner); - } - - return result; -} - -template -py::list pyList1DFromRatCellWrap2D(const T& values) -{ - py::list result; - - for (int32_T idx0{0}; idx0 < values.size(1); idx0++) { - result.append(pyArrayFromRatArray2d(values[idx0].f1)); - } - - return result; -} - -template -py::list pyList1DFromRatCellWrap1D(const T& values) -{ - py::list result; - - for (int32_T idx0{0}; idx0 < values.size(0); idx0++) { - result.append(pyArrayFromRatArray2d(values[idx0].f1)); - } - - return result; -} - -template -py::list pyList2dFromRatCellWrap(const T& values) -{ - py::list result; - int32_T idx {0}; - for (int32_T idx0{0}; idx0 < values.size(0); idx0++) { - py::list inner; - for (int32_T idx1{0}; idx1 < values.size(1); idx1++) { - idx = idx0 + values.size(0) * idx1; - inner.append(pyArrayFromRatArray2d(values[idx].f1)); - } - result.append(inner); - } - - return result; -} - -template -py::list pyListFromBoundedCellWrap(const T& values) -{ - py::list result; - - for (int32_T idx0{0}; idx0 < values.size(1); idx0++) { - auto array = py::array_t({values[idx0].f1.size[0]}); - std::memcpy(array.request().ptr, values[idx0].f1.data, array.nbytes()); - - result.append(array); - } - - return result; -} - -template -py::array_t pyArray1dFromBoundedArray(const T& array) -{ - auto result_array = py::array_t({array.size[0]}); - std::memcpy(result_array.request().ptr, array.data, result_array.nbytes()); - - return result_array; -} - -template -py::array_t pyArray2dFromBoundedArray(const T& array) -{ - auto result_array = py::array_t({array.size[0], array.size[1]}); - std::memcpy(result_array.request().ptr, array.data, result_array.nbytes()); - - return result_array; -} - -py::array_t pyArrayFromRatArray3d(coder::array array) -{ - auto result_array = py::array_t({array.size(0), array.size(1), array.size(2)}); - std::memcpy(result_array.request().ptr, array.data(), result_array.nbytes()); - - return result_array; -} - OutputResult OutputResultFromStruct5T(const RAT::struct5_T result) { // Copy problem to output @@ -1309,6 +591,27 @@ BayesResults bayesResultsFromStruct8T(const RAT::struct8_T results) return bayesResults; } +const std::string docsRATMain = R"(Entry point for the main reflectivity computation. + +Parameters +---------- +problem_def : Rat.rat_core.ProblemDefinition + The project input for the RAT calculation. +limits : RATapi.rat_core.Limits + Min and max values for each parameter defined in the problem definition. +control : RATapi.rat_core.Control + The controls object for the RAT calculation. + +Returns +------- +out_problem_def : Rat.rat_core.ProblemDefinition + The project input with the updated fit values. +results : Rat.rat_core.OutputResult + The results from a RAT calculation. +bayes_result : Rat.rat_core.BayesResults + The extra results if RAT calculation is Bayesian. +)"; + py::tuple RATMain(const ProblemDefinition& problem_def, const Limits& limits, const Control& control) { RAT::struct0_T problem_def_struct = createStruct0(problem_def); @@ -1327,21 +630,42 @@ py::tuple RATMain(const ProblemDefinition& problem_def, const Limits& limits, co bayesResultsFromStruct8T(bayesResults)); } +const std::string docsMakeSLDProfileXY = R"(Creates the profiles for the SLD plots + +Parameters +---------- +bulk_in : float + Bulk in value for contrast. +bulk_out : float + Bulk out value for contrast. +ssub : float + Substrate roughness. +layers : np.ndarray[np.float] + Array of parameters for each layer in the contrast. +number_of_repeats : int, default: 1 + Number of times the layers are repeated. + +Returns +------- +sld_profile : np.ndarray[np.float] + Computed SLD profile +)"; + py::array_t makeSLDProfileXY(real_T bulk_in, real_T bulk_out, real_T ssub, const py::array_t &layers, - real_T number_of_layers, - real_T repeats) + int number_of_repeats=DEFAULT_NREPEATS) { coder::array out; coder::array layers_array = pyArrayToRatArray2d(layers); + py::buffer_info buffer_info = layers.request(); RAT::makeSLDProfileXY(bulk_in, bulk_out, ssub, layers_array, - number_of_layers, - repeats, + buffer_info.shape[0], + number_of_repeats, out); return pyArrayFromRatArray2d(out); @@ -1387,13 +711,13 @@ PYBIND11_MODULE(rat_core, m) { .def("invoke", overload_cast_&, std::vector&>()(&DylibEngine::invoke), py::arg("xdata"), py::arg("param")); - py::class_(m, "PredictionIntervals") + py::class_(m, "PredictionIntervals", docsPredictionIntervals.c_str()) .def(py::init<>()) .def_readwrite("reflectivity", &PredictionIntervals::reflectivity) .def_readwrite("sld", &PredictionIntervals::sld) .def_readwrite("sampleChi", &PredictionIntervals::sampleChi); - py::class_(m, "PlotEventData") + py::class_(m, "PlotEventData", docsPlotEventData.c_str()) .def(py::init<>()) .def_readwrite("reflectivity", &PlotEventData::reflectivity) .def_readwrite("shiftedData", &PlotEventData::shiftedData) @@ -1430,7 +754,7 @@ PYBIND11_MODULE(rat_core, m) { return evt; })); - py::class_(m, "ProgressEventData") + py::class_(m, "ProgressEventData", docsProgressEventData.c_str()) .def(py::init<>()) .def_readwrite("message", &ProgressEventData::message) .def_readwrite("percent", &ProgressEventData::percent) @@ -1452,7 +776,7 @@ PYBIND11_MODULE(rat_core, m) { return evt; })); - py::class_(m, "ConfidenceIntervals") + py::class_(m, "ConfidenceIntervals", docsConfidenceIntervals.c_str()) .def(py::init<>()) .def_readwrite("percentile95", &ConfidenceIntervals::percentile95) .def_readwrite("percentile65", &ConfidenceIntervals::percentile65) @@ -1480,14 +804,14 @@ PYBIND11_MODULE(rat_core, m) { .def_readwrite("storeOutput", &DreamParams::storeOutput) .def_readwrite("R", &DreamParams::R); - py::class_(m, "NestedSamplerOutput") + py::class_(m, "NestedSamplerOutput", docsNestedSamplerOutput.c_str()) .def(py::init<>()) .def_readwrite("logZ", &NestedSamplerOutput::logZ) .def_readwrite("logZErr", &NestedSamplerOutput::logZErr) .def_readwrite("nestSamples", &NestedSamplerOutput::nestSamples) .def_readwrite("postSamples", &NestedSamplerOutput::postSamples); - py::class_(m, "DreamOutput") + py::class_(m, "DreamOutput", docsDreamOutput.c_str()) .def(py::init<>()) .def_readwrite("allChains", &DreamOutput::allChains) .def_readwrite("outlierChains", &DreamOutput::outlierChains) @@ -1498,7 +822,7 @@ PYBIND11_MODULE(rat_core, m) { .def_readwrite("R_stat", &DreamOutput::R_stat) .def_readwrite("CR", &DreamOutput::CR); - py::class_(m, "BayesResults") + py::class_(m, "BayesResults", docsBayesResults.c_str()) .def(py::init<>()) .def_readwrite("predictionIntervals", &BayesResults::predictionIntervals) .def_readwrite("confidenceIntervals", &BayesResults::confidenceIntervals) @@ -1507,12 +831,12 @@ PYBIND11_MODULE(rat_core, m) { .def_readwrite("nestedSamplerOutput", &BayesResults::nestedSamplerOutput) .def_readwrite("chain", &BayesResults::chain); - py::class_(m, "Calculation") + py::class_(m, "Calculation", docsCalculation.c_str()) .def(py::init<>()) .def_readwrite("chiValues", &Calculation::chiValues) .def_readwrite("sumChi", &Calculation::sumChi); - py::class_(m, "ContrastParams") + py::class_(m, "ContrastParams", docsContrastParams.c_str()) .def(py::init<>()) .def_readwrite("scalefactors", &ContrastParams::scalefactors) .def_readwrite("bulkIn", &ContrastParams::bulkIn) @@ -1520,7 +844,7 @@ PYBIND11_MODULE(rat_core, m) { .def_readwrite("subRoughs", &ContrastParams::subRoughs) .def_readwrite("resample", &ContrastParams::resample); - py::class_(m, "OutputResult") + py::class_(m, "OutputResult", docsOutputResult.c_str()) .def(py::init<>()) .def_readwrite("reflectivity", &OutputResult::reflectivity) .def_readwrite("simulation", &OutputResult::simulation) @@ -1535,7 +859,7 @@ PYBIND11_MODULE(rat_core, m) { .def_readwrite("fitParams", &OutputResult::fitParams) .def_readwrite("fitNames", &OutputResult::fitNames); - py::class_(m, "NameStore") + py::class_(m, "NameStore", docsNameStore.c_str()) .def(py::init<>()) .def_readwrite("params", &NameStore::params) .def_readwrite("backgroundParams", &NameStore::backgroundParams) @@ -1572,7 +896,7 @@ PYBIND11_MODULE(rat_core, m) { return names; })); - py::class_(m, "Checks") + py::class_(m, "Checks", docsChecks.c_str()) .def(py::init<>()) .def_readwrite("params", &Checks::params) .def_readwrite("backgroundParams", &Checks::backgroundParams) @@ -1607,7 +931,7 @@ PYBIND11_MODULE(rat_core, m) { return chk; })); - py::class_(m, "Limits") + py::class_(m, "Limits", docsLimits.c_str()) .def(py::init<>()) .def_readwrite("params", &Limits::params) .def_readwrite("backgroundParams", &Limits::backgroundParams) @@ -1642,7 +966,7 @@ PYBIND11_MODULE(rat_core, m) { return lim; })); - py::class_(m, "Control") + py::class_(m, "Control", docsControl.c_str()) .def(py::init<>()) .def_readwrite("parallel", &Control::parallel) .def_readwrite("procedure", &Control::procedure) @@ -1723,7 +1047,7 @@ PYBIND11_MODULE(rat_core, m) { return ctrl; })); - py::class_(m, "ProblemDefinition") + py::class_(m, "ProblemDefinition", docsProblemDefinition.c_str()) .def(py::init<>()) .def_readwrite("TF", &ProblemDefinition::TF) .def_readwrite("resample", &ProblemDefinition::resample) @@ -1859,7 +1183,8 @@ PYBIND11_MODULE(rat_core, m) { return p; })); - m.def("RATMain", &RATMain, "Entry point for the main reflectivity computation."); + m.def("RATMain", &RATMain, docsRATMain.c_str(), py::arg("problem_def"), py::arg("limits"), py::arg("control")); - m.def("makeSLDProfileXY", &makeSLDProfileXY, "Creates the profiles for the SLD plots"); + m.def("makeSLDProfileXY", &makeSLDProfileXY, docsMakeSLDProfileXY.c_str(), + py::arg("bulk_in"), py::arg("bulk_out"), py::arg("ssub"), py::arg("layers"), py::arg("number_of_repeats") = DEFAULT_NREPEATS); } diff --git a/setup.py b/setup.py index a21faae8..84e7c9ce 100644 --- a/setup.py +++ b/setup.py @@ -26,6 +26,7 @@ pybind11.get_include(), pybind11.get_include(True), "cpp/RAT/", + "cpp/includes/", ], language="c++", ), diff --git a/tests/test_plotting.py b/tests/test_plotting.py index a1954834..78ddd08c 100644 --- a/tests/test_plotting.py +++ b/tests/test_plotting.py @@ -73,8 +73,8 @@ def bayes_fig(request) -> plt.figure: @pytest.mark.parametrize("fig", [False, True], indirect=True) -def test_figure_axis_formating(fig: plt.figure) -> None: - """Tests the axis formating of the figure.""" +def test_figure_axis_formatting(fig: plt.figure) -> None: + """Tests the axis formatting of the figure.""" ref_plot = fig.axes[0] sld_plot = fig.axes[1] @@ -151,20 +151,17 @@ def test_sld_profile_function_call(mock: MagicMock) -> None: assert mock.call_args_list[0].args[0] == 2.07e-06 assert mock.call_args_list[0].args[1] == 6.28e-06 assert mock.call_args_list[0].args[2] == 0.0 - assert mock.call_args_list[0].args[4] == 82 - assert mock.call_args_list[0].args[5] == 1.0 + assert mock.call_args_list[0].args[4] == 1 assert mock.call_args_list[1].args[0] == 2.07e-06 assert mock.call_args_list[1].args[1] == 1.83e-06 assert mock.call_args_list[1].args[2] == 0.0 - assert mock.call_args_list[1].args[4] == 128 - assert mock.call_args_list[1].args[5] == 1.0 + assert mock.call_args_list[1].args[4] == 1 assert mock.call_args_list[2].args[0] == 2.07e-06 assert mock.call_args_list[2].args[1] == -5.87e-07 assert mock.call_args_list[2].args[2] == 0.0 - assert mock.call_args_list[2].args[4] == 153 - assert mock.call_args_list[2].args[5] == 1.0 + assert mock.call_args_list[2].args[4] == 1 @patch("RATapi.utils.plotting.makeSLDProfileXY") @@ -181,20 +178,17 @@ def test_live_plot(mock: MagicMock) -> None: assert mock.call_args_list[0].args[0] == 2.07e-06 assert mock.call_args_list[0].args[1] == 6.28e-06 assert mock.call_args_list[0].args[2] == 0.0 - assert mock.call_args_list[0].args[4] == 82 - assert mock.call_args_list[0].args[5] == 1.0 + assert mock.call_args_list[0].args[4] == 1 assert mock.call_args_list[1].args[0] == 2.07e-06 assert mock.call_args_list[1].args[1] == 1.83e-06 assert mock.call_args_list[1].args[2] == 0.0 - assert mock.call_args_list[1].args[4] == 128 - assert mock.call_args_list[1].args[5] == 1.0 + assert mock.call_args_list[1].args[4] == 1 assert mock.call_args_list[2].args[0] == 2.07e-06 assert mock.call_args_list[2].args[1] == -5.87e-07 assert mock.call_args_list[2].args[2] == 0.0 - assert mock.call_args_list[2].args[4] == 153 - assert mock.call_args_list[2].args[5] == 1.0 + assert mock.call_args_list[2].args[4] == 1 @patch("RATapi.utils.plotting.plot_ref_sld_helper") diff --git a/tests/test_wrappers.py b/tests/test_wrappers.py index a50127e4..2ced5ad0 100644 --- a/tests/test_wrappers.py +++ b/tests/test_wrappers.py @@ -27,7 +27,7 @@ def test_matlab_wrapper() -> None: mocked_engine.demo.return_value = ([2], 5) result = handle([1], [2], [3], 0) assert result == ([2], 5) - assert wrapper.engine.demo.call_args[0] == ([1], [2], [3], 1, -1) + assert wrapper.engine.demo.call_args[0] == ([1], [2], [3], 1) mocked_engine.demo.assert_called_once() mocked_engine.demo.return_value = ([3, 1], 7)