diff --git a/CMakeLists.txt b/CMakeLists.txt index 7441196..6402e60 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -45,6 +45,7 @@ option(MaCh3Tutorial_PROFILING_ENABLED "Will run with gperftools" OFF) option(Use_External_MaCh3 "Use an externally installed MaCh3 instead of CPM" OFF) option(MaCh3Tutorial_DEBUG_ENABLED "Use debug options" OFF) option(MaCh3Tutorial_LOW_MEMORY_STRUCTS_ENABLED "Use floats" OFF) +option(PYTHON_ENABLED "Build python package" OFF) # KS: Here we try to find tag matching tutorial version. If we can't find one then use develop # This will allow to grab tutorial for tagged MaCh3 version without a need of manually changing version @@ -113,6 +114,7 @@ if(NOT MaCh3_FOUND) "Prob3ppLinear_ENABLED ON" "NuFastLinear_ENABLED ON" "CUDAProb3_ENABLED ON" + "MaCh3_PYTHON_ENABLED ${PYTHON_ENABLED}" ) # MaCh3 sets up cmessage for us, so can't cmessage until we've set up MaCh3 cmessage(STATUS "Didn't find MaCh3, will checkout the ${MaCh3_Branch} branch from mach3-software/MaCh3") @@ -146,6 +148,51 @@ add_subdirectory(Tutorial) add_subdirectory(SplinesTutorial) add_subdirectory(SamplesTutorial) +# MOVE THIS BLOCK HERE - BEFORE add_subdirectory(python) +# After MaCh3 is added/found, determine MaCh3_PREFIX +if(NOT DEFINED MaCh3_PREFIX OR MaCh3_PREFIX STREQUAL "") + # Try to detect if MaCh3 was fetched by CPM + if(EXISTS "${CMAKE_BINARY_DIR}/_deps/mach3-src") + set(MaCh3_PREFIX "${CMAKE_BINARY_DIR}/_deps/mach3-src") + cmessage(STATUS "Detected CPM-fetched MaCh3: ${MaCh3_PREFIX}") + # Try to get from MaCh3 package config + elseif(DEFINED MaCh3_DIR) + get_filename_component(MaCh3_PREFIX "${MaCh3_DIR}" DIRECTORY) + get_filename_component(MaCh3_PREFIX "${MaCh3_PREFIX}" DIRECTORY) + cmessage(STATUS "Detected installed MaCh3 from MaCh3_DIR: ${MaCh3_PREFIX}") + # Try to get from target properties + elseif(TARGET MaCh3::All) + get_target_property(MACH3_INCLUDE_DIRS MaCh3::All INTERFACE_INCLUDE_DIRECTORIES) + if(MACH3_INCLUDE_DIRS) + list(GET MACH3_INCLUDE_DIRS 0 FIRST_INCLUDE_DIR) + get_filename_component(MaCh3_PREFIX "${FIRST_INCLUDE_DIR}" DIRECTORY) + cmessage(STATUS "Detected MaCh3 from target include dirs: ${MaCh3_PREFIX}") + endif() + endif() + + # Final validation + if(NOT DEFINED MaCh3_PREFIX OR MaCh3_PREFIX STREQUAL "") + cmessage(FATAL_ERROR "Could not determine MaCh3_PREFIX. Please set it manually with -DMaCh3_PREFIX=") + endif() + + if(NOT EXISTS "${MaCh3_PREFIX}") + cmessage(FATAL_ERROR "MaCh3_PREFIX does not exist: ${MaCh3_PREFIX}") + endif() + + # Cache it for use in subdirectories + set(MaCh3_PREFIX "${MaCh3_PREFIX}" CACHE PATH "MaCh3 installation or source directory" FORCE) + cmessage(STATUS "Set MaCh3_PREFIX: ${MaCh3_PREFIX}") +endif() + +# NOW add the python subdirectory - MaCh3_PREFIX is already set +if(PYTHON_ENABLED) + cmessage(STATUS "PYTHON IS ENABLED!!") + set(CMAKE_INSTALL_RPATH "$ORIGIN/../lib") + cmessage(STATUS "CMAKE_INSTALL_RPATH is ${CMAKE_INSTALL_RPATH}") + set(PYMACH3_PATH ${CMAKE_BINARY_DIR}/python/pyMaCh3) + add_subdirectory(python) +endif() + configure_file(cmake/Templates/setup.MaCh3Tutorial.sh.in "${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/setup.MaCh3Tutorial.sh" @ONLY) install(FILES "${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/setup.MaCh3Tutorial.sh" DESTINATION ${CMAKE_INSTALL_PREFIX}/bin) diff --git a/SamplesTutorial/CMakeLists.txt b/SamplesTutorial/CMakeLists.txt index f5e5d21..0af636a 100755 --- a/SamplesTutorial/CMakeLists.txt +++ b/SamplesTutorial/CMakeLists.txt @@ -1,6 +1,6 @@ set(HEADERS - SampleHandlerTutorial.h - StructsTutorial.h + ${CMAKE_CURRENT_SOURCE_DIR}/SampleHandlerTutorial.h + ${CMAKE_CURRENT_SOURCE_DIR}/StructsTutorial.h ) add_library(SamplesTutorial SHARED diff --git a/SplinesTutorial/CMakeLists.txt b/SplinesTutorial/CMakeLists.txt index c2b6c81..a3e1f51 100644 --- a/SplinesTutorial/CMakeLists.txt +++ b/SplinesTutorial/CMakeLists.txt @@ -1,4 +1,4 @@ -set(HEADERS BinnedSplinesTutorial.h +set(HEADERS ${CMAKE_CURRENT_SOURCE_DIR}/BinnedSplinesTutorial.h ) add_library(BinnedSplinesTutorial SHARED diff --git a/cmake/Templates/setup.MaCh3Tutorial.sh.in b/cmake/Templates/setup.MaCh3Tutorial.sh.in index db250e6..bebcf4f 100755 --- a/cmake/Templates/setup.MaCh3Tutorial.sh.in +++ b/cmake/Templates/setup.MaCh3Tutorial.sh.in @@ -108,3 +108,28 @@ else fi unset SETUPDIR + +# Python bindings configuration +PYTHON_ENABLED=@MaCh3_PYTHON_ENABLED@ + +if [ "${PYTHON_ENABLED}" = "ON" ] || [ "${PYTHON_ENABLED}" = "TRUE" ] || [ "${PYTHON_ENABLED}" = "1" ]; then + # Python bindings are enabled + PYMACH3_PATH="@PYMACH3_PATH@" + + if [ -d "${PYMACH3_PATH}" ]; then + if [ -z "${PYTHONPATH}" ]; then + export PYTHONPATH="${PYMACH3_PATH}" + else + # Prepend to existing PYTHONPATH, avoiding duplicates + if [[ ":${PYTHONPATH}:" != *":${PYMACH3_PATH}:"* ]]; then + export PYTHONPATH="${PYMACH3_PATH}:${PYTHONPATH}" + fi + fi + echo "✓ Added pyMaCh3 to PYTHONPATH: ${PYMACH3_PATH}" + else + echo "⚠️ Warning: pyMaCh3 directory not found at ${PYMACH3_PATH}" + echo " Python bindings may not have been built successfully" + fi +else + echo "Python bindings not enabled in this build" +fi \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..3955591 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,20 @@ +[build-system] +requires = ["scikit-build-core", "pybind11"] +build-backend = "scikit_build_core.build" + +[project] +name = "pyMaCh3_tutorial" +version = "0.0.1" +authors = [{ name = "The MaCh3 Collaboration" }] +dependencies = [ + "numpy", + "matplotlib", + "pandas", + "pyMC", +] + +[tool.scikit-build] +wheel.install-dir = "pyMaCh3_tutorial" + +[tool.scikit-build.cmake] +args = ["-DPYTHON_ENABLED=ON"] diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt new file mode 100644 index 0000000..0459d38 --- /dev/null +++ b/python/CMakeLists.txt @@ -0,0 +1,10 @@ + +include(${MaCh3_PREFIX}/cmake/Modules/pyMaCh3.cmake) + +setup_pyMaCh3( + INSTALL_DIR "NONE" + BINDING_FILES pyMaCh3.cpp + LINK_TARGETS + BinnedSplinesTutorial + SamplesTutorial +) diff --git a/python/PyMaCh3Example.ipynb b/python/PyMaCh3Example.ipynb new file mode 100644 index 0000000..7162d93 --- /dev/null +++ b/python/PyMaCh3Example.ipynb @@ -0,0 +1,653 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "e072aa46", + "metadata": {}, + "source": [ + "# Welcome to pyMaCh3!\n", + "\n", + "This notebook is meant to provide a brief introduction to the pyMaCh3 library. It will demonstrate basic functionality and examples of how to use the bindings that have been made from the C++ code of MaCh3. As more bindings are added then more examples can be added.\n", + "\n", + "The first step is setting up your environment to find and use pyMaCh3. For python to find the library that has been made, you need to setup your PYTHONPATH environment variable to include the directory where the library is installed. If you are using a conda environment then I would recommend adding that PYTHONPATH to your envrionment to make it easier to find the library." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "413d0743", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "['BinnedSplinesTutorial', 'SamplesTutorial']\n", + "MaCh3_ROOT env var: /home/ewan/software/MaCh3/MaCh3-core-python/build\n", + "Initializing SampleHandlerTutorial bindings... \n" + ] + } + ], + "source": [ + "import os\n", + "\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "\n", + "# Now import pyMaCh3\n", + "import pyMaCh3_tutorial as m3\n", + "\n", + "print(f'MaCh3_ROOT env var: {os.environ[\"MaCh3_ROOT\"]}')\n" + ] + }, + { + "cell_type": "markdown", + "id": "80b5935b", + "metadata": {}, + "source": [ + "# First steps\n", + "\n", + "Now that pyMaCh3 is imported we can do some basic operations on MaCh3 objects. We can try creating parameter and sample handler objects and make some nice plots from them." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4db4aded", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[2026-03-02 15:01:24.555] [info] [ParameterHandlerBase.cpp:30] Constructing instance of ParameterHandler using\n", + "[2026-03-02 15:01:24.555] [info] [ParameterHandlerBase.cpp:34] ../TutorialConfigs/CovObjs/SystematicModel.yaml\n", + "[2026-03-02 15:01:24.555] [info] [ParameterHandlerBase.cpp:36] as an input\n", + "[2026-03-02 15:01:24.555] [info] [ParameterHandlerBase.cpp:39] Principal component analysis but given the threshold for the principal components to be less than 0, or greater than (or equal to) 1. This will not work\n", + "[2026-03-02 15:01:24.555] [info] [ParameterHandlerBase.cpp:40] Please specify a number between 0 and 1\n", + "[2026-03-02 15:01:24.555] [info] [ParameterHandlerBase.cpp:41] You specified: \n", + "[2026-03-02 15:01:24.555] [info] [ParameterHandlerBase.cpp:42] Am instead calling the usual non-PCA constructor...\n", + "[2026-03-02 15:01:24.689] [info] [ParameterTunes.cpp:50] Found 2 tunes:\n", + "[2026-03-02 15:01:24.689] [info] [ParameterTunes.cpp:52] Tune 0 Generated\n", + "[2026-03-02 15:01:24.689] [info] [ParameterTunes.cpp:52] Tune 1 PostND\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerBase.cpp:371] Created covariance matrix from files: \n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerBase.cpp:373] ../TutorialConfigs/CovObjs/SystematicModel.yaml \n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerBase.cpp:375] ----------------\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerBase.cpp:376] Found 10 systematics parameters in total\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerBase.cpp:377] ----------------\n", + "[2026-03-02 15:01:24.689] [warning] [ParameterHandlerGeneric.cpp:281] Spline knot capping enabled with bounds [0.5, 2]. For reliable fits, consider modifying the input generation instead.\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:454] #################################################\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:455] Printing ParameterHandlerGeneric:\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:478] ============================================================================================================================================================\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:479] # | Name | Prior | Error | Lower | Upper | StepScale | SampleNames | Type \n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:480] ------------------------------------------------------------------------------------------------------------------------------------------------------------\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:490] 0 | Norm_Param_0 | 1 | +/- 0.11 | 0 | 999 | 0.2 | Tutorial_*, ND2137 | Norm \n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:490] 1 | Norm_Param_1 | 1 | +/- 0.18 | 0 | 999 | 0.2 | Tutorial_*, ND2137 | Norm \n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:490] 2 | Norm_Param_2 | 1 | +/- 0.40 | 0 | 9999 | 0.2 | Tutorial_*, ND2137 | Norm \n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:490] 3 | BinnedSplineParam1 | 1 | +/- 1.00 | 0 | 2 | 0.2 | Tutorial_*, ND2137 | Spline \n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:490] 4 | BinnedSplineParam2 | 1 | +/- 1.00 | 0 | 2 | 0.2 | Tutorial_*, ND2137 | Spline \n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:490] 5 | BinnedSplineParam3 | 1 | +/- 1.00 | 0 | 2 | 0.2 | Tutorial_*, ND2137 | Spline \n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:490] 6 | BinnedSplineParam4 | 1 | +/- 1.00 | 0 | 2 | 0.2 | Tutorial_*, ND2137 | Spline \n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:490] 7 | BinnedSplineParam5 | 1 | +/- 1.00 | 0 | 2 | 0.2 | Tutorial_*, ND2137 | Spline \n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:490] 8 | EResLep | 0 | +/- 0.01 | -1 | 1 | 0.2 | Tutorial_* | Functional\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:490] 9 | EResTot | 0 | +/- 0.01 | -1 | 1 | 0.2 | Tutorial_* | Functional\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:492] ============================================================================================================================================================\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:499] Normalisation parameters: 3\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:505] ┌────┬──────────┬────────────────────────────────────────┬────────────────────┬────────────────────┬────────────────────┐\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:506] │# │Global # │Name │Int. mode │Target │pdg │\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:507] ├────┼──────────┼────────────────────────────────────────┼────────────────────┼───���────────────────┼────────────────────┤\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:532] │0 │0 │Norm_Param_0 │1 │12 16 │all │\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:532] │1 │1 │Norm_Param_1 │1 │12 16 │all │\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:532] │2 │2 │Norm_Param_2 │1 │12 16 │all │\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:536] └────┴──────────┴────────────────────────────────────────┴────────────────────┴────────────────────┴────────────────────┘\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:539] Normalisation parameters KinematicCuts information\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:540] ┌────┬──────────┬────────────────────────────────────────┬────────────────────┬────────────────────────────────────────┐\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:541] │# │Global # │Name │KinematicCut │Value │\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:542] ├────┼──────────┼────────────────────────────────────────┼────────────────────┼────────────────────────────────────────┤\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:557] │0 │0 │Norm_Param_0 │TrueQ2 │0.25 0.50 │\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:557] │1 │1 │Norm_Param_1 │TrueQ2 │0.50 1.00 │\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:557] │2 │2 │Norm_Param_2 │TrueQ2 │1.00 9999999.90 │\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:562] ��────┴──────────┴────────────────────────────────────────┴────────────────────┴────────────────────────────────────────┘\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:571] Spline parameters: 5\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:573] =====================================================================================================================================================================\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:574] # | Name | Spline Name | Spline Interpolation | Low Knot Bound | Up Knot Bound | \n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:575] ---------------------------------------------------------------------------------------------------------------------------------------------------------------------\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:580] 0 | BinnedSplineParam1 | mysyst1 | TSpline3 | 0.5 | 2 | \n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:580] 1 | BinnedSplineParam2 | mysyst2 | Linear | -9999 | 9999 | \n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:580] 2 | BinnedSplineParam3 | mysyst3 | Monotonic | -9999 | 9999 | \n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:580] 3 | BinnedSplineParam4 | mysyst4 | Akima | -9999 | 9999 | \n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:580] 4 | BinnedSplineParam5 | mysyst5 | Linear | -9999 | 9999 | \n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:587] =====================================================================================================================================================================\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:593] Functional parameters: 2\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:595] ┌────┬──────────┬────────────────────────────────────────┐\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:596] │# │Global # │Name │\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:597] ├────┼──────────┼────────────────────────────────────────┤\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:601] │0 │8 │EResLep │\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:601] │1 │9 │EResTot │\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:603] └────┴──────────┴────────────────────────────────────────┘\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:609] Oscillation parameters: 0\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:633] Printing parameter groups\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:636] Found 10: Xsec params\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:469] Finished\n", + "[2026-03-02 15:01:24.689] [info] [ParameterHandlerGeneric.cpp:470] #################################################\n" + ] + }, + { + "ename": "", + "evalue": "", + "output_type": "error", + "traceback": [ + "\u001b[1;31mThe Kernel crashed while executing code in the current cell or a previous cell. \n", + "\u001b[1;31mPlease review the code in the cell(s) to identify a possible cause of the failure. \n", + "\u001b[1;31mClick here for more info. \n", + "\u001b[1;31mView Jupyter log for further details." + ] + } + ], + "source": [ + "# Create a parameter handler object from the SystematicModel.yaml file\n", + "parameter_handler = m3.parameters.ParameterHandlerGeneric([\"../TutorialConfigs/CovObjs/SystematicModel.yaml\"])" + ] + }, + { + "cell_type": "markdown", + "id": "848aaa2b", + "metadata": {}, + "source": [ + "Voila! We now have a parameterHandler object to play with in python. For fun, let's throw from this matrix, save the parameter values in a dataframe and then plot them using matplotlib." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b73be13a", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Number of parameters: 10\n" + ] + } + ], + "source": [ + "# Check how many parameters there are in the parameter handler\n", + "num_parameters = parameter_handler.get_n_pars()\n", + "print(f\"Number of parameters: {num_parameters}\")\n", + "\n", + "n_toys = 1000\n", + "\n", + "# Create an array of the parameter names\n", + "parameter_names = [f\"{parameter_handler.get_fancy_par_name(i)}\" for i in range(num_parameters)]\n", + "\n", + "# Create a dataframe to store parameter values with the parameter names as the columns and the number of toys as the index\n", + "parameters_df = pd.DataFrame(columns=parameter_names)\n", + "\n", + "for toy in range(n_toys):\n", + " # throw the parameter_handler object\n", + " parameter_handler.propose_step()\n", + " parameters_df.loc[toy] = parameter_handler.get_proposal_array()\n" + ] + }, + { + "cell_type": "markdown", + "id": "5910c56d", + "metadata": {}, + "source": [ + "Now we can simply use some of the functionality from pandas to draw a scatter plot of the parameters." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d76d6c65", + "metadata": {}, + "outputs": [], + "source": [ + "# Plot a scatter matrix of the parameters\n", + "Axes = pd.plotting.scatter_matrix(parameters_df, alpha=0.2, figsize=(24,12), diagonal='kde')\n", + "\n", + "# This just makes the y-axis labels somewheat readable (or at least not terribly overlapping)\n", + "[plt.setp(item.yaxis.get_label(), 'size', 6) for item in Axes.ravel()]\n", + "\n", + "# Let's show the plot!\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "dcae3eac", + "metadata": {}, + "source": [ + "As you can see, this all looks perfectly sensible and you can even see the inbuild prior correlations between some parameters.\n", + "\n", + "We can now move on to creating a sample handler object." + ] + }, + { + "cell_type": "markdown", + "id": "e5a25dc9", + "metadata": {}, + "source": [ + "# Sample Handler example\n", + "\n", + "We will create a sample handler object and produce some spectra from a few toys we threw." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "43c76314", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Help on module pyMaCh3_tutorial._pyMaCh3.samples in pyMaCh3_tutorial._pyMaCh3:\n", + "\n", + "NAME\n", + " pyMaCh3_tutorial._pyMaCh3.samples - This is a Python binding of MaCh3s C++ based samples library.\n", + "\n", + "CLASSES\n", + " pybind11_builtins.pybind11_object(builtins.object)\n", + " SampleHandlerBase\n", + " SampleHandlerFD\n", + " SampleHandlerTutorial(SampleHandlerFD, SampleHandlerBase)\n", + " TestStatistic\n", + " \n", + " class SampleHandlerBase(pybind11_builtins.pybind11_object)\n", + " | Method resolution order:\n", + " | SampleHandlerBase\n", + " | pybind11_builtins.pybind11_object\n", + " | builtins.object\n", + " | \n", + " | Methods defined here:\n", + " | \n", + " | __init__(...)\n", + " | __init__(self: pyMaCh3_tutorial._pyMaCh3.samples.SampleHandlerBase) -> None\n", + " | \n", + " | get_bin_LLH(...)\n", + " | get_bin_LLH(self: pyMaCh3_tutorial._pyMaCh3.samples.SampleHandlerBase, data: float, mc: float, w2: float) -> float\n", + " | \n", + " | Get the LLH for a bin by comparing the data and MC. The result depends on having previously set the test statistic using :py:meth:`pyMaCh3.samples.SampleHandlerBase.set_test_stat` \n", + " | :param data: The data content of the bin. \n", + " | :param mc: The mc content of the bin \n", + " | :param w2: The Sum(w_{i}^2) (sum of weights squared) in the bin, which is sigma^2_{MC stats}\n", + " | \n", + " | get_likelihood(...)\n", + " | get_likelihood(self: pyMaCh3_tutorial._pyMaCh3.samples.SampleHandlerBase) -> float\n", + " | \n", + " | Get the sample likelihood at the current point in your model space. You will need to override this.\n", + " | \n", + " | reweight(...)\n", + " | reweight(self: pyMaCh3_tutorial._pyMaCh3.samples.SampleHandlerBase) -> None\n", + " | \n", + " | reweight the MC events in this sample. You will need to override this.\n", + " | \n", + " | set_test_stat(...)\n", + " | set_test_stat(self: pyMaCh3_tutorial._pyMaCh3.samples.SampleHandlerBase, test_stat: pyMaCh3_tutorial._pyMaCh3.samples.TestStatistic) -> None\n", + " | \n", + " | Set the test statistic that should be used when calculating likelihoods. \n", + " | :param test_stat: The new test statistic to use\n", + " | \n", + " | ----------------------------------------------------------------------\n", + " | Static methods inherited from pybind11_builtins.pybind11_object:\n", + " | \n", + " | __new__(*args, **kwargs) from pybind11_builtins.pybind11_type\n", + " | Create and return a new object. See help(type) for accurate signature.\n", + " \n", + " class SampleHandlerFD(SampleHandlerBase)\n", + " | Method resolution order:\n", + " | SampleHandlerFD\n", + " | SampleHandlerBase\n", + " | pybind11_builtins.pybind11_object\n", + " | builtins.object\n", + " | \n", + " | Methods defined here:\n", + " | \n", + " | __init__(...)\n", + " | __init__(self: pyMaCh3_tutorial._pyMaCh3.samples.SampleHandlerFD, mc_version: str, xsec_cov: ParameterHandlerGeneric) -> None\n", + " | \n", + " | This should never be called directly as SampleHandlerFD is an abstract base class. \n", + " | However when creating a derived class, in the __init__() method, you should call the parent constructor i.e. this one by doing:: \n", + " | \n", + " | super(, self).__init__(*args)\n", + " | \n", + " | get_data_hist(...)\n", + " | get_data_hist(self: pyMaCh3_tutorial._pyMaCh3.samples.SampleHandlerFD, Dimension: int) -> tuple\n", + " | \n", + " | Get Data histogram as numpy arrays.\n", + " | For 1D: Returns (contents, edges)\n", + " | For 2D: Returns (contents, edgesX, edgesY)\n", + " | where contents is shape (nbinsY, nbinsX) for 2D\n", + " | \n", + " | get_mc_hist(...)\n", + " | get_mc_hist(self: pyMaCh3_tutorial._pyMaCh3.samples.SampleHandlerFD, sample: int) -> tuple\n", + " | \n", + " | Get MC histogram as numpy arrays.\n", + " | For 1D: Returns (contents, edges)\n", + " | For 2D: Returns (contents, edgesX, edgesY)\n", + " | where contents is shape (nbinsY, nbinsX) for 2D\n", + " | \n", + " | get_w2_hist(...)\n", + " | get_w2_hist(self: pyMaCh3_tutorial._pyMaCh3.samples.SampleHandlerFD, Dimension: int) -> tuple\n", + " | \n", + " | Get W2 histogram as numpy arrays.\n", + " | For 1D: Returns (contents, edges)\n", + " | For 2D: Returns (contents, edgesX, edgesY)\n", + " | where contents is shape (nbinsY, nbinsX) for 2D\n", + " | \n", + " | ----------------------------------------------------------------------\n", + " | Methods inherited from SampleHandlerBase:\n", + " | \n", + " | get_bin_LLH(...)\n", + " | get_bin_LLH(self: pyMaCh3_tutorial._pyMaCh3.samples.SampleHandlerBase, data: float, mc: float, w2: float) -> float\n", + " | \n", + " | Get the LLH for a bin by comparing the data and MC. The result depends on having previously set the test statistic using :py:meth:`pyMaCh3.samples.SampleHandlerBase.set_test_stat` \n", + " | :param data: The data content of the bin. \n", + " | :param mc: The mc content of the bin \n", + " | :param w2: The Sum(w_{i}^2) (sum of weights squared) in the bin, which is sigma^2_{MC stats}\n", + " | \n", + " | get_likelihood(...)\n", + " | get_likelihood(self: pyMaCh3_tutorial._pyMaCh3.samples.SampleHandlerBase) -> float\n", + " | \n", + " | Get the sample likelihood at the current point in your model space. You will need to override this.\n", + " | \n", + " | reweight(...)\n", + " | reweight(self: pyMaCh3_tutorial._pyMaCh3.samples.SampleHandlerBase) -> None\n", + " | \n", + " | reweight the MC events in this sample. You will need to override this.\n", + " | \n", + " | set_test_stat(...)\n", + " | set_test_stat(self: pyMaCh3_tutorial._pyMaCh3.samples.SampleHandlerBase, test_stat: pyMaCh3_tutorial._pyMaCh3.samples.TestStatistic) -> None\n", + " | \n", + " | Set the test statistic that should be used when calculating likelihoods. \n", + " | :param test_stat: The new test statistic to use\n", + " | \n", + " | ----------------------------------------------------------------------\n", + " | Static methods inherited from pybind11_builtins.pybind11_object:\n", + " | \n", + " | __new__(*args, **kwargs) from pybind11_builtins.pybind11_type\n", + " | Create and return a new object. See help(type) for accurate signature.\n", + " \n", + " class SampleHandlerTutorial(SampleHandlerFD, SampleHandlerBase)\n", + " | Method resolution order:\n", + " | SampleHandlerTutorial\n", + " | SampleHandlerFD\n", + " | SampleHandlerBase\n", + " | pybind11_builtins.pybind11_object\n", + " | builtins.object\n", + " | \n", + " | Methods defined here:\n", + " | \n", + " | __init__(...)\n", + " | __init__(*args, **kwargs)\n", + " | Overloaded function.\n", + " | \n", + " | 1. __init__(self: pyMaCh3_tutorial._pyMaCh3.samples.SampleHandlerTutorial, mc_version: str, xsec_cov: ParameterHandlerGeneric) -> None\n", + " | \n", + " | Create SampleHandlerTutorial without oscillation handler\n", + " | \n", + " | 2. __init__(self: pyMaCh3_tutorial._pyMaCh3.samples.SampleHandlerTutorial, mc_version: str, xsec_cov: ParameterHandlerGeneric, osc_cov: OscillationHandler = None) -> None\n", + " | \n", + " | Create SampleHandlerTutorial with oscillation handler\n", + " | \n", + " | ----------------------------------------------------------------------\n", + " | Methods inherited from SampleHandlerFD:\n", + " | \n", + " | get_data_hist(...)\n", + " | get_data_hist(self: pyMaCh3_tutorial._pyMaCh3.samples.SampleHandlerFD, Dimension: int) -> tuple\n", + " | \n", + " | Get Data histogram as numpy arrays.\n", + " | For 1D: Returns (contents, edges)\n", + " | For 2D: Returns (contents, edgesX, edgesY)\n", + " | where contents is shape (nbinsY, nbinsX) for 2D\n", + " | \n", + " | get_mc_hist(...)\n", + " | get_mc_hist(self: pyMaCh3_tutorial._pyMaCh3.samples.SampleHandlerFD, sample: int) -> tuple\n", + " | \n", + " | Get MC histogram as numpy arrays.\n", + " | For 1D: Returns (contents, edges)\n", + " | For 2D: Returns (contents, edgesX, edgesY)\n", + " | where contents is shape (nbinsY, nbinsX) for 2D\n", + " | \n", + " | get_w2_hist(...)\n", + " | get_w2_hist(self: pyMaCh3_tutorial._pyMaCh3.samples.SampleHandlerFD, Dimension: int) -> tuple\n", + " | \n", + " | Get W2 histogram as numpy arrays.\n", + " | For 1D: Returns (contents, edges)\n", + " | For 2D: Returns (contents, edgesX, edgesY)\n", + " | where contents is shape (nbinsY, nbinsX) for 2D\n", + " | \n", + " | ----------------------------------------------------------------------\n", + " | Methods inherited from SampleHandlerBase:\n", + " | \n", + " | get_bin_LLH(...)\n", + " | get_bin_LLH(self: pyMaCh3_tutorial._pyMaCh3.samples.SampleHandlerBase, data: float, mc: float, w2: float) -> float\n", + " | \n", + " | Get the LLH for a bin by comparing the data and MC. The result depends on having previously set the test statistic using :py:meth:`pyMaCh3.samples.SampleHandlerBase.set_test_stat` \n", + " | :param data: The data content of the bin. \n", + " | :param mc: The mc content of the bin \n", + " | :param w2: The Sum(w_{i}^2) (sum of weights squared) in the bin, which is sigma^2_{MC stats}\n", + " | \n", + " | get_likelihood(...)\n", + " | get_likelihood(self: pyMaCh3_tutorial._pyMaCh3.samples.SampleHandlerBase) -> float\n", + " | \n", + " | Get the sample likelihood at the current point in your model space. You will need to override this.\n", + " | \n", + " | reweight(...)\n", + " | reweight(self: pyMaCh3_tutorial._pyMaCh3.samples.SampleHandlerBase) -> None\n", + " | \n", + " | reweight the MC events in this sample. You will need to override this.\n", + " | \n", + " | set_test_stat(...)\n", + " | set_test_stat(self: pyMaCh3_tutorial._pyMaCh3.samples.SampleHandlerBase, test_stat: pyMaCh3_tutorial._pyMaCh3.samples.TestStatistic) -> None\n", + " | \n", + " | Set the test statistic that should be used when calculating likelihoods. \n", + " | :param test_stat: The new test statistic to use\n", + " | \n", + " | ----------------------------------------------------------------------\n", + " | Static methods inherited from pybind11_builtins.pybind11_object:\n", + " | \n", + " | __new__(*args, **kwargs) from pybind11_builtins.pybind11_type\n", + " | Create and return a new object. See help(type) for accurate signature.\n", + " \n", + " class TestStatistic(pybind11_builtins.pybind11_object)\n", + " | Method resolution order:\n", + " | TestStatistic\n", + " | pybind11_builtins.pybind11_object\n", + " | builtins.object\n", + " | \n", + " | Methods defined here:\n", + " | \n", + " | __eq__(...)\n", + " | __eq__(self: object, other: object) -> bool\n", + " | \n", + " | __getstate__(...)\n", + " | __getstate__(self: object) -> int\n", + " | \n", + " | __hash__(...)\n", + " | __hash__(self: object) -> int\n", + " | \n", + " | __index__(...)\n", + " | __index__(self: pyMaCh3_tutorial._pyMaCh3.samples.TestStatistic) -> int\n", + " | \n", + " | __init__(...)\n", + " | __init__(self: pyMaCh3_tutorial._pyMaCh3.samples.TestStatistic, value: int) -> None\n", + " | \n", + " | __int__(...)\n", + " | __int__(self: pyMaCh3_tutorial._pyMaCh3.samples.TestStatistic) -> int\n", + " | \n", + " | __ne__(...)\n", + " | __ne__(self: object, other: object) -> bool\n", + " | \n", + " | __repr__(...)\n", + " | __repr__(self: object) -> str\n", + " | \n", + " | __setstate__(...)\n", + " | __setstate__(self: pyMaCh3_tutorial._pyMaCh3.samples.TestStatistic, state: int) -> None\n", + " | \n", + " | __str__(...)\n", + " | __str__(self: object) -> str\n", + " | \n", + " | ----------------------------------------------------------------------\n", + " | Readonly properties defined here:\n", + " | \n", + " | __members__\n", + " | \n", + " | name\n", + " | name(self: object) -> str\n", + " | \n", + " | value\n", + " | \n", + " | ----------------------------------------------------------------------\n", + " | Data and other attributes defined here:\n", + " | \n", + " | Barlow_Beeston = \n", + " | \n", + " | Dembinski_Abdelmottele = \n", + " | \n", + " | Ice_Cube = \n", + " | \n", + " | N_Test_Statistics = \n", + " | \n", + " | Pearson = \n", + " | \n", + " | Poisson = \n", + " | \n", + " | ----------------------------------------------------------------------\n", + " | Static methods inherited from pybind11_builtins.pybind11_object:\n", + " | \n", + " | __new__(*args, **kwargs) from pybind11_builtins.pybind11_type\n", + " | Create and return a new object. See help(type) for accurate signature.\n", + "\n", + "FILE\n", + " (built-in)\n", + "\n", + "\n" + ] + } + ], + "source": [ + "#help(m3.samples)\n", + "\n", + "sample_handler = m3.samples.SampleHandlerTutorial(str(\"../TutorialConfigs/Samples/SampleHandler_Tutorial.yaml\"), parameter_handler)" + ] + }, + { + "cell_type": "markdown", + "id": "ba663db3", + "metadata": {}, + "source": [ + "As you can see we've now create a sample handler object and you can see the usual C++ stdout from MaCh3. We can now produce a histogram of the MC prediction from MaCh3." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b3335d2e", + "metadata": {}, + "outputs": [ + { + "ename": "NameError", + "evalue": "name 'sample_handler' is not defined", + "output_type": "error", + "traceback": [ + "\u001b[31m---------------------------------------------------------------------------\u001b[39m", + "\u001b[31mNameError\u001b[39m Traceback (most recent call last)", + "\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[6]\u001b[39m\u001b[32m, line 1\u001b[39m\n\u001b[32m----> \u001b[39m\u001b[32m1\u001b[39m \u001b[43msample_handler\u001b[49m.reweight()\n\u001b[32m 2\u001b[39m mc_prediction, bin_edges = sample_handler.get_mc_hist(\u001b[32m1\u001b[39m)\n\u001b[32m 4\u001b[39m \u001b[38;5;66;03m# count the number of non-zero entries in the MC prediction\u001b[39;00m\n", + "\u001b[31mNameError\u001b[39m: name 'sample_handler' is not defined" + ] + } + ], + "source": [ + "sample_handler.reweight()\n", + "mc_prediction, bin_edges = sample_handler.get_mc_hist(1)\n", + "\n", + "# count the number of non-zero entries in the MC prediction\n", + "non_zero_entries = np.count_nonzero(mc_prediction)\n", + "print(f\"Number of non-zero entries in MC prediction: {non_zero_entries}\")\n", + "\n", + "#Now pass this to matplotlib to draw a histogram\n", + "plt.figure(figsize=(10, 6))\n", + "plt.hist(bin_edges[:-1], bins=bin_edges, weights=mc_prediction, histtype='step', \n", + " linewidth=2, label=\"MC prediction\")\n", + "plt.xlabel(\"Energy [GeV]\", fontsize=12)\n", + "plt.ylabel(\"Events\", fontsize=12)\n", + "plt.legend(fontsize=11)\n", + "plt.grid(True, alpha=0.3)\n", + "plt.title(\"MaCh3 Tutorial - MC Prediction\", fontsize=14)\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "d851d335", + "metadata": {}, + "source": [ + "Finally, let's bring the two bits of code together and produce many spectra from the sample handler from the parameter handler throws we previously made." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b72eb7ae", + "metadata": {}, + "outputs": [], + "source": [ + "# Produce many spectra using the parameter values in the parameters_df\n", + "mc_prediction_integral = np.zeros(n_toys)\n", + "\n", + "for toy in range(n_toys):\n", + " # Set the parameter values from the dataframe\n", + " parameter_handler.set_parameters(parameters_df.iloc[toy].values)\n", + " # Re-weight the sample handler\n", + " sample_handler.reweight()\n", + " # Get the MC prediction and integrate it\n", + " mc_prediction, _ = sample_handler.get_mc_hist(1)\n", + " mc_prediction_integral[toy] = np.sum(mc_prediction)\n", + "\n", + "# Plot the integrals of the MC predictions\n", + "plt.figure(figsize=(10, 6))\n", + "plt.hist(mc_prediction_integral, bins=20, histtype='step', \n", + " linewidth=2, label=\"MC prediction integral\")\n", + "plt.xlabel(\"MC prediction integral\", fontsize=12)\n", + "plt.ylabel(\"Events\", fontsize=12)\n", + "plt.legend(fontsize=11)" + ] + }, + { + "cell_type": "markdown", + "id": "488c5bda", + "metadata": {}, + "source": [ + "And just like that we have a prior uncertainty on the total number of events in our prediction. This is a very simple example but it might inspire more complicated python scripts." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "env", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.0" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/python/pyMaCh3.cpp b/python/pyMaCh3.cpp new file mode 100644 index 0000000..488ddc6 --- /dev/null +++ b/python/pyMaCh3.cpp @@ -0,0 +1,45 @@ +#include "python/pyMaCh3.h" + +#include "SamplesTutorial/SampleHandlerTutorial.h" +#include + +namespace py = pybind11; + +class MaCh3TutorialPyBinder : public MaCh3PyBinder { + + public: + + void initSamplesExperiment(py::module &m_samples){ + + std::cout << "Initializing SampleHandlerTutorial bindings... " << std::endl; + + // Now add SampleHandlerTutorial to the same submodule + py::class_(m_samples, "SampleHandlerTutorial") + // Constructor with 2 arguments (no oscillation handler) + .def(py::init([](const std::string& mc_version, ParameterHandlerGeneric* xsec_cov) { + return new SampleHandlerTutorial(mc_version, xsec_cov, nullptr); + }), + "Create SampleHandlerTutorial without oscillation handler", + py::arg("mc_version"), + py::arg("xsec_cov") + ) + // Constructor with 3 arguments (with oscillation handler) + .def(py::init([](const std::string& mc_version, + ParameterHandlerGeneric* xsec_cov, + OscillationHandler* osc_cov) { + std::shared_ptr osc_ptr; + if (osc_cov != nullptr) { + osc_ptr = std::shared_ptr(osc_cov, [](OscillationHandler*){}); + } + return new SampleHandlerTutorial(mc_version, xsec_cov, osc_ptr); + }), + "Create SampleHandlerTutorial with oscillation handler", + py::arg("mc_version"), + py::arg("xsec_cov"), + py::arg("osc_cov") = nullptr + ) + ; + } +}; + +MAKE_PYMACH3_MDULE( MaCh3TutorialPyBinder )