diff --git a/.github/workflows/build_wheels.yml b/.github/workflows/build_wheels.yml
index c3a53fd1d..695bce351 100644
--- a/.github/workflows/build_wheels.yml
+++ b/.github/workflows/build_wheels.yml
@@ -26,41 +26,77 @@ jobs:
     steps:
       - uses: actions/checkout@v4
 
-      # Install Fortran compiler based on OS
-      - name: Install dependencies
-        run: |
-          if [ "${{ matrix.plat.os }}" = "ubuntu-latest" ]; then
-            sudo apt-get update
-            sudo apt-get install -y gfortran
-          else
-            brew update
-            brew install gfortran
-          fi
+      - name: Restore compiler cache
+        uses: actions/cache@v4
+        with:
+          path: _build/ccache
+          key: pyalps-ccache-${{ runner.os }}-${{ matrix.plat.arch }}-${{ hashFiles('src/**', 'bindings/python/**', 'applications/**', 'tool/maxent*') }}
+          restore-keys: pyalps-ccache-${{ runner.os }}-${{ matrix.plat.arch }}-
 
       - name: Build wheels
-        uses: pypa/cibuildwheel@v2.22.0
+        uses: pypa/cibuildwheel@v3.4.1
+        with:
+          package-dir: bindings/python/pyalps
         env:
-          CIBW_BUILD: cp39-* cp310-* cp311-* cp312-* cp313-*
+          CIBW_BUILD: cp310-* cp311-* cp312-* cp313-* cp314-*
           CIBW_ARCHS: ${{ matrix.plat.arch }}
           CIBW_ARCHS_MACOS: ${{ matrix.plat.arch }}
-
-          # Set Fortran compiler for all platforms
-          CIBW_ENVIRONMENT: "FC=gfortran"
-
-          # macOS-specific settings
+          CIBW_ENVIRONMENT: >
+            ALPS_DIR=$(pwd)/_build/cibw-install/share/alps
+            CCACHE_DIR=$(pwd)/_build/ccache
+            CCACHE_NAMESPACE=pyalps-wheel
+            CMAKE_ARGS="-DCMAKE_CXX_COMPILER_LAUNCHER=ccache"
+          # manylinux (AlmaLinux/glibc) uses dnf and glibc's SunRPC/XDR + execinfo.
+          # musllinux (Alpine/musl) has neither: install libtirpc for the system
+          # XDR path (ALPS_HAVE_RPC_XDR_H) and disable the execinfo backtrace.
+          CIBW_BEFORE_ALL_LINUX: >
+            if command -v dnf >/dev/null 2>&1;
+            then dnf install -y ccache cmake hdf5-devel lapack-devel ninja-build; EXTRA=;
+            else apk add --no-cache ccache ninja-is-really-ninja hdf5-dev lapack-dev libtirpc-dev;
+            ln -sf /usr/include/tirpc/rpc /usr/include/rpc;
+            ln -sf /usr/include/tirpc/netconfig.h /usr/include/netconfig.h;
+            export CXXFLAGS="-DALPS_NGS_NO_STACKTRACE";
+            EXTRA="-DALPS_HAVE_RPC_XDR_H=1 -DCMAKE_CXX_STANDARD_LIBRARIES=-ltirpc";
+            fi &&
+            cmake -S {project} -B {project}/_build/cibw-alps -G Ninja
+            -DCMAKE_BUILD_TYPE=Release
+            -DCMAKE_INSTALL_PREFIX={project}/_build/cibw-install
+            -DCMAKE_CXX_COMPILER_LAUNCHER=ccache
+            -DALPS_BUILD_LIBS_ONLY=ON
+            -DALPS_BUILD_TESTS=OFF
+            -DALPS_BUILD_EXAMPLES=OFF
+            -DALPS_BUILD_APPLICATIONS=OFF
+            -DALPS_ENABLE_MPI=OFF
+            $EXTRA &&
+            cmake --build {project}/_build/cibw-alps --target install -j2
+          CIBW_BEFORE_ALL_MACOS: >
+            brew install ccache cmake hdf5 ninja &&
+            cmake -S {project} -B {project}/_build/cibw-alps -G Ninja
+            -DCMAKE_BUILD_TYPE=Release
+            -DCMAKE_INSTALL_PREFIX={project}/_build/cibw-install
+            -DCMAKE_CXX_COMPILER_LAUNCHER=ccache
+            -DALPS_BUILD_LIBS_ONLY=ON
+            -DALPS_BUILD_TESTS=OFF
+            -DALPS_BUILD_EXAMPLES=OFF
+            -DALPS_BUILD_APPLICATIONS=OFF
+            -DALPS_ENABLE_MPI=OFF
+            -DHDF5_ROOT=${{ matrix.plat.homebrew }}/opt/hdf5 &&
+            cmake --build {project}/_build/cibw-alps --target install -j2
           CIBW_ENVIRONMENT_MACOS: >
-            Boost_ROOT_DIR=/Users/runner/work/ALPS/ALPS/boost_1_87_0
+            ALPS_DIR=$(pwd)/_build/cibw-install/share/alps
+            CCACHE_DIR=$(pwd)/_build/ccache
+            CCACHE_NAMESPACE=pyalps-wheel
+            CMAKE_ARGS="-DCMAKE_CXX_COMPILER_LAUNCHER=ccache"
             MACOSX_DEPLOYMENT_TARGET=${{ matrix.plat.target }}
             CXXFLAGS="-stdlib=libc++"
-#          CIBW_ENVIRONMENT: >
-
-        # env:
-        #   CIBW_SOME_OPTION: value
-        #    ...
-        # with:
-        #   package-dir: .
-        #   output-dir: wheelhouse
-        #   config-file: "{package}/pyproject.toml"
+          CIBW_REPAIR_WHEEL_COMMAND_LINUX: >
+            auditwheel repair -w {dest_dir} {wheel} &&
+            auditwheel show {dest_dir}/*.whl
+          CIBW_REPAIR_WHEEL_COMMAND_MACOS: >
+            delocate-wheel --require-archs {delocate_archs} -w {dest_dir} -v {wheel} &&
+            delocate-listdeps --all {dest_dir}/*.whl
+          CIBW_TEST_REQUIRES: pytest
+          CIBW_TEST_COMMAND: pytest -q {project}/test/pyalps
 
       - uses: actions/upload-artifact@v4
         with:
@@ -75,7 +111,7 @@ jobs:
       - uses: actions/checkout@v4
 
       - name: Build sdist
-        run: pipx run build --sdist
+        run: pipx run build --sdist --outdir dist bindings/python/pyalps
 
       - uses: actions/upload-artifact@v4
         with:
diff --git a/CMakeLists.txt b/CMakeLists.txt
index cbcb17674..ad7203a89 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -28,12 +28,6 @@ if (NOT_BUILD_SHARED_LIBS)
 else (NOT_BUILD_SHARED_LIBS)
 	option(BUILD_SHARED_LIBS "Build shared libraries" ON)
 endif (NOT_BUILD_SHARED_LIBS)
-if (NOT_ALPS_BUILD_PYTHON)
-	option(ALPS_BUILD_PYTHON "Build ALPS python extentions" OFF)
-else (NOT_ALPS_BUILD_PYTHON)
-	option(ALPS_BUILD_PYTHON "Build ALPS python extentions" ON)
-endif (NOT_ALPS_BUILD_PYTHON)
-
 option(ALPS_INSTALL_DEPRECATED "Install deprecated headers for ALPS" OFF)
 option(ALPS_BUILD_DEVELOPER_TOOLS "Build tools used by developers to maintain ALPS" OFF)
 option(ALPS_ENABLE_OPENMP "Enable OpenMP parallelization" OFF)
@@ -47,7 +41,6 @@ option(ALPS_LINK_BOOST_TEST "Link Boost.test pre-built library" OFF)
 option(ALPS_INSTALL_BOOST_TEST "Install Boost Test framework" OFF)
 option(ALPS_NGS_USE_NEW_ALEA "Use the ALPS ngs-alea instead of just alea" OFF)
 option(ALPS_NGS_OPENMPI_ULFM "Userlevel failure mitigation is available" OFF)
-option(ALPS_PYTHON_WHEEL "Build python wheel" OFF)
 option(ALPS_BUILD_LIBS_ONLY "Build only libraries" OFF)
 
 
@@ -57,7 +50,6 @@ mark_as_advanced(ALPS_LINK_BOOST_TEST)
 mark_as_advanced(ALPS_INSTALL_BOOST_TEST)
 mark_as_advanced(ALPS_NGS_USE_NEW_ALEA)
 mark_as_advanced(ALPS_NGS_OPENMPI_ULFM)
-mark_as_advanced(ALPS_PYTHON_WHEEL)
 mark_as_advanced(ALPS_BUILD_LIBS_ONLY)
 
 option(ALPS_USE_MKL_PARALLEL "Use parallel version of MKL" OFF)
@@ -66,7 +58,6 @@ mark_as_advanced(ALPS_USE_MKL_PARALLEL)
 SET (APPLICATIONS_CAN_BE_BUILT ON)
 
 set(ALPS_BOOST_LIBRARY_NAME "boost" CACHE STRING "name of the boost library")
-set(ALPS_BOOST_PYTHON_LIBRARY_NAME "boost_python" CACHE STRING "name of the boost library")
 set(ALPS_ENABLE_MPI ON CACHE BOOL "Enable MPI Parallelization")
 set(ALPS_INSTALL_HEADERS  ON CACHE BOOL "Install headers for ALPS and all dependent libraries")
 set(ALPS_BUILD_EXAMPLES ON CACHE BOOL "Build ALPS examples")
@@ -75,11 +66,6 @@ set(ALPS_BUILD_APPLICATIONS ${APPLICATIONS_CAN_BE_BUILT} CACHE BOOL "Build ALPS
 
 mark_as_advanced(ALPS_BOOST_LIBRARY_NAME)
 
-if(ALPS_PYTHON_WHEEL)
-  set(ALPS_BUILD_APPLICATIONS ON)
-  set(ALPS_ENABLE_MPI OFF)
-endif()
-
 if(ALPS_BUILD_LIBS_ONLY)
   set(ALPS_BUILD_APPLICATIONS OFF)
 endif()
@@ -129,10 +115,6 @@ list(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/config)
 # set default CMAKE_INSTALL_PREFIX
 ######################################################################
 
-if (ALPS_BUILD_PYTHON)
-  find_package(PythonMod REQUIRED) # COMPONENTS Interpreter Development.Module )
-endif (ALPS_BUILD_PYTHON)
-
 ######################################################################
 # Version information
 ######################################################################
@@ -276,24 +258,8 @@ IF(HDF5_IS_PARALLEL)
   INCLUDE_DIRECTORIES(${MPI_INCLUDE_PATH})
 ENDIF(HDF5_IS_PARALLEL)
 
-# python
-
-
-if (ALPS_BUILD_PYTHON)
-  set(PYTHON_SCRIPTDIR "${CMAKE_INSTALL_PREFIX}/lib/python")
-  find_package(PythonMod)
-endif (ALPS_BUILD_PYTHON)
-
-IF (PYTHONLIBS_FOUND AND ALPS_BUILD_PYTHON)
-  include_directories(${PYTHON_NUMPY_INCLUDE_DIR})
-    MESSAGE (STATUS "Numpy include in ${PYTHON_NUMPY_INCLUDE_DIR}")
-  SET(ALPS_HAVE_PYTHON ON)
-  INCLUDE_DIRECTORIES(${PYTHON_INCLUDE_DIRS})
-  set(BUILD_BOOST_PYTHON TRUE)
-ELSE (PYTHONLIBS_FOUND AND ALPS_BUILD_PYTHON)
-  set(BUILD_BOOST_PYTHON OFF)
-  SET(ALPS_HAVE_PYTHON OFF)
-ENDIF (PYTHONLIBS_FOUND AND ALPS_BUILD_PYTHON)
+# Python bindings are built by the standalone scikit-build-core project in
+# bindings/python/pyalps.  The C++ SDK deliberately has no Python dependency.
 
 # Boost Libraries
 find_package(BoostForALPS REQUIRED)
@@ -396,12 +362,6 @@ endif(MPI_INCLUDE_DIR)
 if(HDF5_INCLUDE_DIR)
   list(APPEND ALPS_EXTRA_INCLUDE_DIRS ${HDF5_INCLUDE_DIR})
 endif(HDF5_INCLUDE_DIR)
-if(PYTHON_INCLUDE_DIRS)
-  list(APPEND ALPS_EXTRA_INCLUDE_DIRS ${PYTHON_INCLUDE_DIRS})
-endif(PYTHON_INCLUDE_DIRS)
-if(PYTHON_NUMPY_INCLUDE_DIR)
-  list(APPEND ALPS_EXTRA_INCLUDE_DIRS ${PYTHON_NUMPY_INCLUDE_DIR})
-endif(PYTHON_NUMPY_INCLUDE_DIR)
 if(Boost_INCLUDE_DIR_CONFIG)
   list(APPEND ALPS_EXTRA_INCLUDE_DIRS ${Boost_INCLUDE_DIR_CONFIG})
 endif(Boost_INCLUDE_DIR_CONFIG)
@@ -433,10 +393,6 @@ endif(LAPACK_LIBRARIES)
 if(HDF5_LIBRARIES)
   list(APPEND ALPS_EXTRA_LIBRARIES ${HDF5_LIBRARIES})
 endif(HDF5_LIBRARIES)
-if(PYTHON_LIBRARY)
-  list(APPEND ALPS_EXTRA_LIBRARIES ${PYTHON_LIBRARY})
-endif(PYTHON_LIBRARY)
-
 configure_file(config/ALPSConfig.cmake.in ${PROJECT_BINARY_DIR}/config/ALPSConfig.cmake @ONLY)
 configure_file(config/include.mk.in ${PROJECT_BINARY_DIR}/config/include.mk)
 
@@ -444,22 +400,17 @@ configure_file(config/include.mk.in ${PROJECT_BINARY_DIR}/config/include.mk)
 # installation
 ######################################################################
 include(InstallRequiredSystemLibraries)
-if(ALPS_PYTHON_WHEEL)
-    install(DIRECTORY lib/pyalps DESTINATION .  COMPONENT python
-          FILES_MATCHING PATTERN "*.py"
-         )
-    install(DIRECTORY lib/xml DESTINATION pyalps COMPONENT xml
-          FILES_MATCHING PATTERN "*.xsl"
-         )
-    install(DIRECTORY ${PROJECT_BINARY_DIR}/lib/xml DESTINATION pyalps COMPONENT xml)
-elseif(ALPS_BUILD_LIBS_ONLY)
-else()
-
 if (ALPS_INSTALL_HEADERS)
   set(ALPS_HEADER_DIR "include")
   install(DIRECTORY src/alps src/boost src/ietl src/mocasito COMPONENT headers DESTINATION ${ALPS_HEADER_DIR}
           FILES_MATCHING PATTERN "*.h" PATTERN "*.hpp" PATTERN "*.ipp"
          )
+  # ALPS' public headers include the vendored Boost.Numeric bindings.  Install
+  # them as part of the C++ SDK so downstream consumers do not need the ALPS
+  # source tree on their include path.
+  install(DIRECTORY bindings/boost COMPONENT headers DESTINATION ${ALPS_HEADER_DIR}
+          FILES_MATCHING PATTERN "*.h" PATTERN "*.hpp" PATTERN "*.ipp"
+         )
   install(FILES ${PROJECT_BINARY_DIR}/src/alps/config.h ${PROJECT_BINARY_DIR}/src/alps/version.h DESTINATION ${ALPS_HEADER_DIR}/alps COMPONENT headers)
   if (ALPS_INSTALL_DEPRECATED)
     install(DIRECTORY legacy/alps COMPONENT deprecated DESTINATION ${ALPS_HEADER_DIR}
@@ -472,12 +423,6 @@ install(DIRECTORY lib/xml DESTINATION ${ALPS_XML_PATH} COMPONENT xml
         FILES_MATCHING PATTERN "*.xsl"
        )
 
-if(ALPS_BUILD_PYTHON AND ALPS_PYTHON_LIB_DEST_ROOT)
-  install(DIRECTORY lib/pyalps DESTINATION ${ALPS_PYTHON_LIB_DEST_ROOT}  COMPONENT python
-        FILES_MATCHING PATTERN "*.py" PATTERN "*.pyc"
-       )
-endif()
-
 install(DIRECTORY ${PROJECT_BINARY_DIR}/lib/xml DESTINATION ${ALPS_XML_PATH} COMPONENT xml)
 
 add_subdirectory(config)
@@ -492,43 +437,24 @@ install(FILES config/UseALPS.cmake
 
 install(FILES LICENSE.txt README.md
         DESTINATION share/alps COMPONENT libraries)
-if(ALPS_BUILD_PYTHON AND ALPS_PYTHON_LIB_DEST_ROOT)
-  string(CONFIGURE [[
-   set(PROJECT_SOURCE_DIR "@PROJECT_SOURCE_DIR@")
-   set(ALPS_PYTHON_LIB_DEST_ROOT "@ALPS_PYTHON_LIB_DEST_ROOT@")
-   message(STATUS "PROJECT_SOURCE_DIR: ${PROJECT_SOURCE_DIR}")
-   message(STATUS "CMAKE_BINARY_DIR: ${CMAKE_BINARY_DIR}")
-   message(STATUS "CMAKE_SOURCE_DIR: ${CMAKE_SOURCE_DIR}")
-   message(STATUS "CMAKE_INSTALL_PREFIX: ${CMAKE_INSTALL_PREFIX}")
-   message(STATUS "ALPS_PYTHON_LIB_DEST_ROOT: ${ALPS_PYTHON_LIB_DEST_ROOT}")
-
-   configure_file(${PROJECT_SOURCE_DIR}/lib/pyalps/pyalps_config.py.in ${CMAKE_INSTALL_PREFIX}/${ALPS_PYTHON_LIB_DEST_ROOT}/pyalps/pyalps_config.py)
-  ]] install_script @ONLY)
-  install(CODE ${install_script})
-endif()
-
-set(PYTHONEXEC python)
-
-endif()
 ######################################################################
 # libraries
 ######################################################################
 
 add_subdirectory(src/boost)
 add_subdirectory(src/alps)
-add_subdirectory(lib/pyalps)
 
 
 ######################################################################
 # programs
 ######################################################################
 
-if (NOT ALPS_PYTHON_WHEEL AND NOT ALPS_BUILD_LIBS_ONLY)
+if (NOT ALPS_BUILD_LIBS_ONLY)
   add_subdirectory(tool)
   if (ALPS_BUILD_EXAMPLES)
     add_subdirectory(example)
   endif (ALPS_BUILD_EXAMPLES)
-endif (NOT ALPS_PYTHON_WHEEL AND NOT ALPS_BUILD_LIBS_ONLY)
+endif (NOT ALPS_BUILD_LIBS_ONLY)
 
 
 ######################################################################
@@ -539,9 +465,9 @@ if (ALPS_BUILD_APPLICATIONS AND NOT ALPS_BUILD_LIBS_ONLY)
   add_subdirectory(applications)
 endif (ALPS_BUILD_APPLICATIONS AND NOT ALPS_BUILD_LIBS_ONLY)
 
-if (ALPS_INCLUDE_TUTORIALS AND NOT ALPS_PYTHON_WHEEL AND NOT ALPS_BUILD_LIBS_ONLY)
+if (ALPS_INCLUDE_TUTORIALS AND NOT ALPS_BUILD_LIBS_ONLY)
   add_subdirectory(tutorials)
-endif (ALPS_INCLUDE_TUTORIALS AND NOT ALPS_PYTHON_WHEEL AND NOT ALPS_BUILD_LIBS_ONLY)
+endif (ALPS_INCLUDE_TUTORIALS AND NOT ALPS_BUILD_LIBS_ONLY)
 
 ######################################################################
 # developer tools
diff --git a/README-py.md b/README-py.md
deleted file mode 100644
index 4bd2da976..000000000
--- a/README-py.md
+++ /dev/null
@@ -1,50 +0,0 @@
-[![ALPS CI/CD](https://github.com/ALPSim/legacy/actions/workflows/build.yml/badge.svg)](https://github.com/ALPSim/legacy/actions/workflows/build.yml)
-
-## Python Algorithms and Libraries for Physics Simulations
-
-This is python packages for `Algorithms and Libraries for Physics Simulations` project. For more information check [README.txt](https://pypi.org/project/pyalps/2.3.3/README.txt).
-
-### Installation instruction from binaries
-
-1. pyALPS can be installed on most Linux and MacOS mcachines from prebuilt biniaries available on [PyPi](https://pypi.org/project/pyalps).
-pyALPS can be installed using `pip` Python package manager:
-
-```
-pip install pyalps
-```
-
-### Installation instruction from sources
-
-1. Prerequisites
-  - CMake > 3.18
-  - Boost sources >= 1.76
-  - BLAS/LAPACK
-  - HDF5
-  - MPI
-  - Python >= 3.9
-    - Python 3.13 requires Boost version 1.87 or later
-    - Earlier versions maybe also work but unsupported
-  - C++ compiler (build has been tested on GCC 10.5 through 14.2)
-  - GNU Make or Ninja build system
-
-You need to download and unpack boost library:
-```
-wget https://archives.boost.io/release/1.86.0/source/boost_1_86_0.tar.gz
-tar -xzf boost_1_86_0.tar.gz
-```
-Here we download `boost v1.86.0`, we have tested ALPS with versions `1.76.0` and `1.86.0`.
-
-2. Downloading and building sources
-```
-git clone https://github.com/alpsim/ALPS ALPS
-cd ALPS
-Boost_SRC_DIR=`pwd`/../boost_1_86_0 python3 -m build --wheel
-```
-This will download the most recent version of ALPS from the github repository, and build pyALPS python package.
-
-3. Installation
-
-Based on the version of the Python used to build pyALPS, the corresponding Python wheel will be created and stored in `dist` subdirectory. It can be installed using `pip`:
-```
-pip install dist/pyalps-<specs>.whl
-```
diff --git a/applications/diag/fulldiag/CMakeLists.txt b/applications/diag/fulldiag/CMakeLists.txt
index aa48df42e..bd1e3785f 100644
--- a/applications/diag/fulldiag/CMakeLists.txt
+++ b/applications/diag/fulldiag/CMakeLists.txt
@@ -27,18 +27,7 @@ if(LAPACK_FOUND)
   add_executable(fulldiag_evaluate fulldiag_evaluate.C)
   target_link_libraries(fulldiag fulldiag_impl)
   target_link_libraries(fulldiag_evaluate fulldiag_impl)
-  if(ALPS_PYTHON_WHEEL)
-  if(APPLE)
-    set_target_properties(fulldiag PROPERTIES INSTALL_RPATH "@loader_path/../lib" )
-    set_target_properties(fulldiag_evaluate PROPERTIES INSTALL_RPATH "@loader_path/../lib" )
-  else(APPLE)
-    set_target_properties(fulldiag PROPERTIES INSTALL_RPATH "$ORIGIN/../lib" )
-    set_target_properties(fulldiag_evaluate PROPERTIES INSTALL_RPATH "$ORIGIN/../lib" )
-  endif(APPLE)
-  install(TARGETS fulldiag fulldiag_evaluate RUNTIME DESTINATION pyalps/bin COMPONENT applications)
-  else()
   install(TARGETS fulldiag fulldiag_evaluate RUNTIME DESTINATION bin COMPONENT applications)
-  endif()
 else(LAPACK_FOUND)
   message(STATUS "fulldiag will not be built since lapack library has not been found")
 endif(LAPACK_FOUND)
diff --git a/applications/diag/sparsediag/CMakeLists.txt b/applications/diag/sparsediag/CMakeLists.txt
index 71f251adc..d3fd8f28a 100644
--- a/applications/diag/sparsediag/CMakeLists.txt
+++ b/applications/diag/sparsediag/CMakeLists.txt
@@ -22,16 +22,7 @@ if(LAPACK_FOUND)
   set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${LAPACK_LINKER_FLAGS}")
   add_executable(sparsediag sparsediag.C factory.C)
   target_link_libraries(sparsediag alps ${LAPACK_LIBRARY} ${BLAS_LIBRARY})
-  if(ALPS_PYTHON_WHEEL)
-  if(APPLE)
-    set_target_properties(sparsediag PROPERTIES INSTALL_RPATH "@loader_path/../lib" )
-  else(APPLE)
-    set_target_properties(sparsediag PROPERTIES INSTALL_RPATH "$ORIGIN/../lib" )
-  endif(APPLE)
-  install(TARGETS sparsediag RUNTIME DESTINATION pyalps/bin COMPONENT applications)
-  else()
   install(TARGETS sparsediag RUNTIME DESTINATION bin COMPONENT applications)
-  endif()
 else(LAPACK_FOUND)
   message(STATUS "sparsediag will not be built since lapack library is not found")
 endif(LAPACK_FOUND)
diff --git a/applications/dmft/qmc/CMakeLists.txt b/applications/dmft/qmc/CMakeLists.txt
index aa95cfe07..394bec0a4 100644
--- a/applications/dmft/qmc/CMakeLists.txt
+++ b/applications/dmft/qmc/CMakeLists.txt
@@ -115,29 +115,11 @@ if(LAPACK_FOUND)
   add_executable(dmft_interaction_expansion_choice dmft_interaction_expansion_choice.C)
   set_property(TARGET dmft_interaction_expansion_choice PROPERTY LABELS dmft)
   add_alps_test(dmft_interaction_expansion_choice)
-if(ALPS_PYTHON_WHEEL)
-  if(APPLE)
-    set_target_properties(dmft PROPERTIES INSTALL_RPATH "@loader_path/../lib" )
-    set_target_properties(hirschfye PROPERTIES INSTALL_RPATH "@loader_path/../lib" )
-    set_target_properties(hybridization PROPERTIES INSTALL_RPATH "@loader_path/../lib" )
-    set_target_properties(interaction PROPERTIES INSTALL_RPATH "@loader_path/../lib" )
-  else(APPLE)
-    set_target_properties(dmft PROPERTIES INSTALL_RPATH "$ORIGIN/../lib" )
-    set_target_properties(hirschfye PROPERTIES INSTALL_RPATH "$ORIGIN/../lib" )
-    set_target_properties(hybridization PROPERTIES INSTALL_RPATH "$ORIGIN/../lib" )
-    set_target_properties(interaction PROPERTIES INSTALL_RPATH "$ORIGIN/../lib" )
-  endif(APPLE)
-  install(TARGETS dmft RUNTIME DESTINATION pyalps/bin COMPONENT applications)
-  install(TARGETS hirschfye  RUNTIME DESTINATION pyalps/bin COMPONENT applications)
-  install(TARGETS hybridization RUNTIME DESTINATION pyalps/bin COMPONENT applications)
-  install(TARGETS interaction RUNTIME DESTINATION pyalps/bin COMPONENT applications)
-else()
   install(TARGETS dmft RUNTIME DESTINATION bin COMPONENT applications)
   install(TARGETS hirschfye  RUNTIME DESTINATION bin COMPONENT applications)
   install(TARGETS hybridization RUNTIME DESTINATION bin COMPONENT applications)
   install(TARGETS interaction RUNTIME DESTINATION bin COMPONENT applications)
   install(FILES hybridization/Documentation/hybdoc.pdf DESTINATION doc)
-endif()
 else(LAPACK_FOUND)
   message(STATUS "dmft will not be built since the lapack library has not been found")
 endif(LAPACK_FOUND)
diff --git a/applications/dmft/qmc/hybridization/Documentation/hybdoc.tex b/applications/dmft/qmc/hybridization/Documentation/hybdoc.tex
index 3ac25a4af..701bc2012 100644
--- a/applications/dmft/qmc/hybridization/Documentation/hybdoc.tex
+++ b/applications/dmft/qmc/hybridization/Documentation/hybdoc.tex
@@ -141,7 +141,7 @@ \subsubsection{Running the solver}
 \subsection{Python interface}
 \label{pythoninterface}
 
-If ALPS is built with Python support (parameter \verb#ALPS_BUILD_PYTHON=ON#), the solver is also built as a Python module. It can directly be called from within a Python script. This provides a flexible framework which allows one to easily set up tasks ranging from calculations for multiple parameters to complex selfconsistency schemes.
+The standalone pyalps wheel includes the solver module when built with application bindings enabled (the default). It can directly be called from within a Python script. This provides a flexible framework which allows one to easily set up tasks ranging from calculations for multiple parameters to complex selfconsistency schemes.
 
 Basic usage of the Python interface is illustrated by the following script, which repeats the previous example for the standalone executable:
 \begin{verbatim}
diff --git a/applications/dmft/qmc/hybridization/hybmain.cpp b/applications/dmft/qmc/hybridization/hybmain.cpp
index df8bd8f85..3a641e0c6 100644
--- a/applications/dmft/qmc/hybridization/hybmain.cpp
+++ b/applications/dmft/qmc/hybridization/hybmain.cpp
@@ -47,11 +47,11 @@ void master_final_tasks(const alps::results_type<hybridization>::type &results,
 int global_mpi_rank;
 
 #ifdef BUILD_PYTHON_MODULE
-//compile it as a python module (requires boost::python library)
-using namespace boost::python;
+#include "dict_to_params.hpp"
+namespace nb = nanobind;
 
-void solve(boost::python::dict parms_){
-  alps::parameters_type<hybridization>::type parms(parms_);
+void solve(nb::dict const & parms_){
+  alps::parameters_type<hybridization>::type parms = pyalps::params_from_dict(parms_);
   std::string output_file = boost::lexical_cast<std::string>(parms["BASENAME"]|"results")+std::string(".out.h5");
 #else
 int main(int argc, char** argv){
@@ -137,12 +137,10 @@ void master_final_tasks(const alps::results_type<hybridization>::type &results,
 }
 
 #ifdef BUILD_PYTHON_MODULE
-BOOST_PYTHON_MODULE(cthyb)
-{
-    def("solve",solve);//define python-callable run method
-};
+NB_MODULE(cthyb, m) {
+    m.def("solve", solve);
+}
 #endif
 
 
 
-
diff --git a/applications/dmft/qmc/interaction_expansion2/main.cpp b/applications/dmft/qmc/interaction_expansion2/main.cpp
index ce7630914..875db3f76 100644
--- a/applications/dmft/qmc/interaction_expansion2/main.cpp
+++ b/applications/dmft/qmc/interaction_expansion2/main.cpp
@@ -44,11 +44,11 @@ void compute_greens_functions(const alps::results_type<HubbardInteractionExpansi
 int global_mpi_rank;
 
 #ifdef BUILD_PYTHON_MODULE
-//compile it as a python module (requires boost::python library)
-using namespace boost::python;
+#include "dict_to_params.hpp"
+namespace nb = nanobind;
 
-void solve(boost::python::dict parms_){
-    alps::parameters_type<HubbardInteractionExpansion>::type parms(parms_);
+void solve(nb::dict const & parms_){
+    alps::parameters_type<HubbardInteractionExpansion>::type parms = pyalps::params_from_dict(parms_);
     std::string output_file = boost::lexical_cast<std::string>(parms["BASENAME"]|"results")+std::string(".out.h5");
 #else
 int main(int argc, char** argv)
@@ -114,9 +114,8 @@ int main(int argc, char** argv)
 }
 
 #ifdef BUILD_PYTHON_MODULE
-    BOOST_PYTHON_MODULE(ctint)
-    {
-        def("solve",solve);//define python-callable run method
-    };
+    NB_MODULE(ctint, m) {
+        m.def("solve", solve);
+    }
 #endif
     
diff --git a/applications/dmrg/dmrg/CMakeLists.txt b/applications/dmrg/dmrg/CMakeLists.txt
index 034689236..18115bcff 100644
--- a/applications/dmrg/dmrg/CMakeLists.txt
+++ b/applications/dmrg/dmrg/CMakeLists.txt
@@ -26,16 +26,7 @@ else(ALPS_LLVM_WORKAROUND)
   set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${LAPACK_LINKER_FLAGS}")
   add_executable(dmrg dmrg.C factory.C)
   target_link_libraries(dmrg alps ${LAPACK_LIBRARY} ${BLAS_LIBRARY})
-if(ALPS_PYTHON_WHEEL)
-  install(TARGETS dmrg RUNTIME DESTINATION pyalps/bin COMPONENT applications)
-  if(APPLE)
-    set_target_properties(dmrg PROPERTIES INSTALL_RPATH "@loader_path/../lib" )
-  else(APPLE)
-    set_target_properties(dmrg PROPERTIES INSTALL_RPATH "$ORIGIN/../lib" )
-  endif(APPLE)
-else()
   install(TARGETS dmrg RUNTIME DESTINATION bin COMPONENT applications)
-endif()
 endif(ALPS_LLVM_WORKAROUND)
 else(LAPACK_FOUND)
   message(STATUS "dmrg will not be built since lapack library is not found")
diff --git a/applications/mc/simple/CMakeLists.txt b/applications/mc/simple/CMakeLists.txt
index 05471ee89..3680f49ae 100644
--- a/applications/mc/simple/CMakeLists.txt
+++ b/applications/mc/simple/CMakeLists.txt
@@ -19,16 +19,7 @@
 
 add_executable(simplemc main.C evaluator.C ising.C xy.C heisenberg.C)
 target_link_libraries(simplemc alps)
-if(ALPS_PYTHON_WHEEL)
-  if(APPLE)
-    set_target_properties(simplemc PROPERTIES INSTALL_RPATH "@loader_path/../lib" )
-  else(APPLE)
-    set_target_properties(simplemc PROPERTIES INSTALL_RPATH "$ORIGIN/../lib" )
-  endif(APPLE)
-  install(TARGETS simplemc RUNTIME DESTINATION pyalps/bin COMPONENT applications)
-else()
-  install(TARGETS simplemc RUNTIME DESTINATION bin COMPONENT applications)
-endif()
+install(TARGETS simplemc RUNTIME DESTINATION bin COMPONENT applications)
 enable_testing()
 add_alps_test(simplemc_ising simplemc ising ising)
 add_alps_test(simplemc_xy simplemc xy xy)
diff --git a/applications/mc/spins/CMakeLists.txt b/applications/mc/spins/CMakeLists.txt
index e89ff4364..dcfb6d807 100644
--- a/applications/mc/spins/CMakeLists.txt
+++ b/applications/mc/spins/CMakeLists.txt
@@ -26,18 +26,7 @@ if(LAPACK_FOUND)
   add_executable(spinmc_evaluate spinmc_evaluate.C)
   target_link_libraries(spinmc spinmc_impl)
   target_link_libraries(spinmc_evaluate spinmc_impl)
-  if(ALPS_PYTHON_WHEEL)
-  if(APPLE)
-    set_target_properties(spinmc PROPERTIES INSTALL_RPATH "@loader_path/../lib" )
-    set_target_properties(spinmc_evaluate PROPERTIES INSTALL_RPATH "@loader_path/../lib" )
-  else(APPLE)
-    set_target_properties(spinmc PROPERTIES INSTALL_RPATH "$ORIGIN/../lib" )
-    set_target_properties(spinmc_evaluate PROPERTIES INSTALL_RPATH "$ORIGIN/../lib" )
-  endif(APPLE)
-  install(TARGETS spinmc spinmc_evaluate RUNTIME DESTINATION pyalps/bin COMPONENT applications)
-  else()
   install(TARGETS spinmc spinmc_evaluate RUNTIME DESTINATION bin COMPONENT applications)
-  endif()
 else(LAPACK_FOUND)
   message(STATUS "spins will not be built since lapack library is not found")
 endif(LAPACK_FOUND)
diff --git a/applications/qmc/checksign/CMakeLists.txt b/applications/qmc/checksign/CMakeLists.txt
index 90511ce0e..7aef774f1 100644
--- a/applications/qmc/checksign/CMakeLists.txt
+++ b/applications/qmc/checksign/CMakeLists.txt
@@ -19,13 +19,4 @@
 
 add_executable(checksign checksign.C)
 target_link_libraries(checksign alps)
-if(ALPS_PYTHON_WHEEL)
-  if(APPLE)
-    set_target_properties(checksign PROPERTIES INSTALL_RPATH "@loader_path/../lib" )
-  else(APPLE)
-    set_target_properties(checksign PROPERTIES INSTALL_RPATH "$ORIGIN/../lib" )
-  endif(APPLE)
-  install(TARGETS checksign RUNTIME DESTINATION pyalps/bin COMPONENT applications)
-else()
-  install(TARGETS checksign RUNTIME DESTINATION bin COMPONENT applications)
-endif()
\ No newline at end of file
+install(TARGETS checksign RUNTIME DESTINATION bin COMPONENT applications)
diff --git a/applications/qmc/dwa/CMakeLists.txt b/applications/qmc/dwa/CMakeLists.txt
index 242f17f49..94c9cc1f5 100644
--- a/applications/qmc/dwa/CMakeLists.txt
+++ b/applications/qmc/dwa/CMakeLists.txt
@@ -22,16 +22,7 @@ if(LAPACK_FOUND)
   set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${LAPACK_LINKER_FLAGS}")
   add_executable(dwa dwa.cpp)
   target_link_libraries(dwa alps ${LAPACK_LIBRARY} ${BLAS_LIBRARY})
-  if(ALPS_PYTHON_WHEEL)
-  if(APPLE)
-    set_target_properties(dwa PROPERTIES INSTALL_RPATH "@loader_path/../lib" )
-  else(APPLE)
-    set_target_properties(dwa PROPERTIES INSTALL_RPATH "$ORIGIN/../lib" )
-  endif(APPLE)
-  install(TARGETS dwa RUNTIME DESTINATION pyalps/bin COMPONENT applications)
-  else()
   install(TARGETS dwa RUNTIME DESTINATION bin COMPONENT applications)
-  endif()
 else(LAPACK_FOUND)
   message(STATUS "dwa will not be built since lapack library has not been found")
 endif(LAPACK_FOUND)
diff --git a/applications/qmc/dwa/bandstructure.hpp b/applications/qmc/dwa/bandstructure.hpp
index b48ce7af4..1071bc5d3 100644
--- a/applications/qmc/dwa/bandstructure.hpp
+++ b/applications/qmc/dwa/bandstructure.hpp
@@ -32,10 +32,6 @@
 #include <vector>
 #include <boost/multi_array.hpp>
 #include <alps/numeric/vector_functions.hpp>
-#include <alps/ngs/boost_python.hpp>
-#include <boost/bind/bind.hpp>
-#include <boost/python/suite/indexing/vector_indexing_suite.hpp>
-#include <alps/python/numpy_array.hpp>
 
 
 // LAPACK Library: dsteqr -- description: diagonalize a real tridiagonal matrix
@@ -47,7 +43,7 @@ class bandstructure
 {
 public:
   bandstructure(double V0_, double lambda_, double a_, double m_, unsigned int L_, int Mmax_=10);
-  bandstructure(boost::python::object V0_, boost::python::object lambda_, double a_, double m_, unsigned int L_, int Mmax_=10); 
+  bandstructure(std::vector<double> const & V0_, std::vector<double> const & lambda_, double a_, double m_, unsigned int L_, int Mmax_=10);
 
   std::vector<double> get_t()    {  if(!is_evaluated)  evaluate();  return t;     }
   double              get_U()    {  if(!is_evaluated)  evaluate();  return U;     }
@@ -132,15 +128,15 @@ bandstructure::bandstructure(double V0_, double lambda_, double a_, double m_, u
   wk2_d = 1.;
 }
 
-bandstructure::bandstructure(boost::python::object V0_, boost::python::object lambda_, double a_, double m_, unsigned int L_, int Mmax_)
+bandstructure::bandstructure(std::vector<double> const & V0_, std::vector<double> const & lambda_, double a_, double m_, unsigned int L_, int Mmax_)
   : is_evaluated (false)
   , L            (L_)
   , Mmax         (Mmax_)
 {
-  alps::python::numpy::convert(V0_, V0);
-  alps::python::numpy::convert(lambda_, lambda);
+  V0     = V0_;
+  lambda = lambda_;
 
-  if (V0.empty() || lambda.empty() || a_ == 0. || m_ == 0. || L == 0) 
+  if (V0.empty() || lambda.empty() || a_ == 0. || m_ == 0. || L == 0)
     boost::throw_exception(std::runtime_error("Illegal initialization parameters for bandstructure class"));
 
   V0.resize    (3, V0.back());
diff --git a/applications/qmc/looper/CMakeLists.txt b/applications/qmc/looper/CMakeLists.txt
index dfdf50c88..1e3ddca04 100644
--- a/applications/qmc/looper/CMakeLists.txt
+++ b/applications/qmc/looper/CMakeLists.txt
@@ -23,16 +23,7 @@ if(LAPACK_FOUND)
   include_directories(${PROJECT_SOURCE_DIR}/applications/qmc/looper)
   add_executable(loop loop.C loop_custom.C loop_model.C path_integral.C sse.C)
   target_link_libraries(loop alps ${LAPACK_LIBRARY} ${BLAS_LIBRARY})
-  if(ALPS_PYTHON_WHEEL)
-  if(APPLE)
-    set_target_properties(loop PROPERTIES INSTALL_RPATH "@loader_path/../lib" )
-  else(APPLE)
-    set_target_properties(loop PROPERTIES INSTALL_RPATH "$ORIGIN/../lib" )
-  endif(APPLE)
-  install(TARGETS loop RUNTIME DESTINATION pyalps/bin COMPONENT applications)
-  else()
   install(TARGETS loop RUNTIME DESTINATION bin COMPONENT applications)
-  endif()
 else(LAPACK_FOUND)
   message(STATUS "loop will not be built since lapack library is not found")
 endif(LAPACK_FOUND)
diff --git a/applications/qmc/qwl/CMakeLists.txt b/applications/qmc/qwl/CMakeLists.txt
index 6474dd5d4..eb7255733 100644
--- a/applications/qmc/qwl/CMakeLists.txt
+++ b/applications/qmc/qwl/CMakeLists.txt
@@ -22,15 +22,4 @@ add_executable(qwl_evaluate qwl_evaluate.C)
 
 target_link_libraries(qwl alps)
 target_link_libraries(qwl_evaluate alps)
-if(ALPS_PYTHON_WHEEL)
-  if(APPLE)
-    set_target_properties(qwl PROPERTIES INSTALL_RPATH "@loader_path/../lib" )
-    set_target_properties(qwl_evaluate PROPERTIES INSTALL_RPATH "@loader_path/../lib" )
-  else(APPLE)
-    set_target_properties(qwl PROPERTIES INSTALL_RPATH "$ORIGIN/../lib" )
-    set_target_properties(qwl_evaluate PROPERTIES INSTALL_RPATH "$ORIGIN/../lib" )
-  endif(APPLE)
-  install(TARGETS qwl qwl_evaluate RUNTIME DESTINATION pyalps/bin COMPONENT applications)
-else()
-  install(TARGETS qwl qwl_evaluate RUNTIME DESTINATION bin COMPONENT applications)
-endif()
+install(TARGETS qwl qwl_evaluate RUNTIME DESTINATION bin COMPONENT applications)
diff --git a/applications/qmc/sse/CMakeLists.txt b/applications/qmc/sse/CMakeLists.txt
index 86e9f2dfc..193621cf3 100644
--- a/applications/qmc/sse/CMakeLists.txt
+++ b/applications/qmc/sse/CMakeLists.txt
@@ -23,29 +23,11 @@ if(LPSolve_FOUND AND NOT MSVC)
                  SSE.Directed.cpp SSE.Initialization.cpp SSE.Measurements.cpp
 		 SSE.Update.cpp SSE.cpp)
   target_link_libraries(dirloop_sse_v1 alps ${LPSolve_LIBRARIES})
-  if(ALPS_PYTHON_WHEEL)
-  if(APPLE)
-    set_target_properties(dirloop_sse_v1 PROPERTIES INSTALL_RPATH "@loader_path/../lib" )
-  else(APPLE)
-    set_target_properties(dirloop_sse_v1 PROPERTIES INSTALL_RPATH "$ORIGIN/../lib" )
-  endif(APPLE)
-  install(TARGETS dirloop_sse_v1 RUNTIME DESTINATION pyalps/bin COMPONENT applications)
-  else()
   install(TARGETS dirloop_sse_v1 RUNTIME DESTINATION bin COMPONENT applications)
-  endif()
 endif(LPSolve_FOUND AND NOT MSVC)
 if (LPSolve_FOUND AND APPLE)
     set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -idirafter /usr/include/malloc")
 endif(LPSolve_FOUND AND APPLE)
 add_executable(dirloop_sse_evaluate evaluate.C)
 target_link_libraries(dirloop_sse_evaluate alps)
-if(ALPS_PYTHON_WHEEL)
-  if(APPLE)
-    set_target_properties(dirloop_sse_evaluate PROPERTIES INSTALL_RPATH "@loader_path/../lib" )
-  else(APPLE)
-    set_target_properties(dirloop_sse_evaluate PROPERTIES INSTALL_RPATH "$ORIGIN/../lib" )
-  endif(APPLE)
-  install(TARGETS dirloop_sse_evaluate RUNTIME DESTINATION pyalps/bin COMPONENT applications)
-else()
-  install(TARGETS dirloop_sse_evaluate RUNTIME DESTINATION bin COMPONENT applications)
-endif()
\ No newline at end of file
+install(TARGETS dirloop_sse_evaluate RUNTIME DESTINATION bin COMPONENT applications)
diff --git a/applications/qmc/sse4/CMakeLists.txt b/applications/qmc/sse4/CMakeLists.txt
index 3ad0ef331..996d86f5e 100644
--- a/applications/qmc/sse4/CMakeLists.txt
+++ b/applications/qmc/sse4/CMakeLists.txt
@@ -19,13 +19,4 @@
 
 add_executable(dirloop_sse main.cc lp_sse.cpp)
 target_link_libraries(dirloop_sse alps)
-if(ALPS_PYTHON_WHEEL)
-  if(APPLE)
-    set_target_properties(dirloop_sse PROPERTIES INSTALL_RPATH "@loader_path/../lib" )
-  else(APPLE)
-    set_target_properties(dirloop_sse PROPERTIES INSTALL_RPATH "$ORIGIN/../lib" )
-  endif(APPLE)
-  install(TARGETS dirloop_sse RUNTIME DESTINATION pyalps/bin COMPONENT applications)
-else()
-  install(TARGETS dirloop_sse RUNTIME DESTINATION bin COMPONENT applications)
-endif()
\ No newline at end of file
+install(TARGETS dirloop_sse RUNTIME DESTINATION bin COMPONENT applications)
diff --git a/applications/qmc/worms/CMakeLists.txt b/applications/qmc/worms/CMakeLists.txt
index 93247388e..823c5fe93 100644
--- a/applications/qmc/worms/CMakeLists.txt
+++ b/applications/qmc/worms/CMakeLists.txt
@@ -25,17 +25,4 @@ add_executable(worm_evaluate evaluate.C)
 
 target_link_libraries(worm worm_impl)
 target_link_libraries(worm_evaluate worm_impl)
-if(ALPS_PYTHON_WHEEL)
-  if(APPLE)
-    set_target_properties(worm PROPERTIES INSTALL_RPATH "@loader_path/../lib" )
-    set_target_properties(worm_evaluate PROPERTIES INSTALL_RPATH "@loader_path/../lib" )
-  else(APPLE)
-    set_target_properties(worm PROPERTIES INSTALL_RPATH "$ORIGIN/../lib" )
-    set_target_properties(worm_evaluate PROPERTIES INSTALL_RPATH "$ORIGIN/../lib" )
-  endif(APPLE)
-  install(TARGETS worm RUNTIME DESTINATION pyalps/bin COMPONENT applications)
-  install(TARGETS worm_evaluate RUNTIME DESTINATION pyalps/bin COMPONENT applications)
-else()
-  install(TARGETS worm RUNTIME DESTINATION bin COMPONENT applications)
-  install(TARGETS worm_evaluate RUNTIME DESTINATION bin COMPONENT applications)
-endif()
\ No newline at end of file
+install(TARGETS worm worm_evaluate RUNTIME DESTINATION bin COMPONENT applications)
diff --git a/bindings/python/pyalps/CMakeLists.txt b/bindings/python/pyalps/CMakeLists.txt
new file mode 100644
index 000000000..ba658bba9
--- /dev/null
+++ b/bindings/python/pyalps/CMakeLists.txt
@@ -0,0 +1,143 @@
+# Copyright (C) 2026 by the ALPS collaboration
+# SPDX-License-Identifier: MIT
+
+cmake_minimum_required(VERSION 3.18)
+project(pyalps LANGUAGES CXX)
+
+option(PYALPS_BUILD_APPLICATIONS "Build optional ALPS application bindings" ON)
+
+set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+set(CMAKE_CXX_EXTENSIONS OFF)
+
+find_package(ALPS REQUIRED CONFIG)
+find_package(Python 3.10 REQUIRED COMPONENTS Interpreter Development.Module)
+
+execute_process(
+  COMMAND "${Python_EXECUTABLE}" -m nanobind --cmake_dir
+  OUTPUT_VARIABLE _nanobind_cmake_dir
+  OUTPUT_STRIP_TRAILING_WHITESPACE
+  COMMAND_ERROR_IS_FATAL ANY)
+list(PREPEND CMAKE_PREFIX_PATH "${_nanobind_cmake_dir}")
+find_package(nanobind 2.10 CONFIG REQUIRED)
+
+get_filename_component(_repo_root "${CMAKE_CURRENT_SOURCE_DIR}/../../.." ABSOLUTE)
+set(_alps_source_root "${_repo_root}")
+if(NOT EXISTS "${_alps_source_root}/tool/maxent.cpp")
+  set(_alps_source_root "${CMAKE_CURRENT_SOURCE_DIR}/_vendor")
+endif()
+set(_bindings "${CMAKE_CURRENT_SOURCE_DIR}/cpp")
+
+link_directories(${ALPS_LIBRARY_DIRS})
+if(ALPS_HDF5_INCLUDE_DIR)
+  get_filename_component(_alps_hdf5_prefix "${ALPS_HDF5_INCLUDE_DIR}" DIRECTORY)
+  link_directories("${_alps_hdf5_prefix}/lib")
+endif()
+
+# Legacy ALPSConfig.cmake serializes the HDF5 imported target name into
+# ALPS_LIBRARIES.  Outside the original build tree that target no longer
+# exists, while the installed library keeps its conventional `hdf5` name.
+set(_pyalps_link_libraries ${ALPS_LIBRARIES})
+list(TRANSFORM _pyalps_link_libraries REPLACE "^hdf5-shared$" "hdf5")
+separate_arguments(_pyalps_link_options NATIVE_COMMAND "${ALPS_EXTRA_LINKER_FLAGS}")
+
+set(_pyalps_targets
+  pyalea_c
+  pymcdata_c
+  pytools_c
+  pyngsparams_c
+  pyngshdf5_c
+  pyngsbase_c
+  pyngsobservable_c
+  pyngsobservables_c
+  pyngsresult_c
+  pyngsresults_c
+  pyngsapi_c
+  pyngsrandom01_c
+  pyngsaccumulator_c)
+
+nanobind_add_module(pyalea_c NB_STATIC "${_bindings}/pyalea.cpp")
+nanobind_add_module(pymcdata_c NB_STATIC "${_bindings}/pymcdata.cpp")
+nanobind_add_module(pytools_c NB_STATIC "${_bindings}/pytools.cpp")
+nanobind_add_module(pyngsparams_c NB_STATIC "${_bindings}/ngs/params.cpp")
+nanobind_add_module(pyngshdf5_c NB_STATIC "${_bindings}/ngs/hdf5.cpp")
+nanobind_add_module(pyngsbase_c NB_STATIC "${_bindings}/ngs/mcbase.cpp")
+nanobind_add_module(pyngsobservable_c NB_STATIC "${_bindings}/ngs/observable.cpp")
+nanobind_add_module(pyngsobservables_c NB_STATIC "${_bindings}/ngs/observables.cpp")
+nanobind_add_module(pyngsresult_c NB_STATIC "${_bindings}/ngs/result.cpp")
+nanobind_add_module(pyngsresults_c NB_STATIC "${_bindings}/ngs/results.cpp")
+nanobind_add_module(pyngsapi_c NB_STATIC "${_bindings}/ngs/api.cpp")
+nanobind_add_module(pyngsrandom01_c NB_STATIC "${_bindings}/ngs/random01.cpp")
+nanobind_add_module(pyngsaccumulator_c NB_STATIC "${_bindings}/ngs/accumulator.cpp")
+
+if(PYALPS_BUILD_APPLICATIONS)
+  if(NOT EXISTS "${_alps_source_root}/tool/maxent.cpp")
+    message(FATAL_ERROR
+      "PYALPS_BUILD_APPLICATIONS requires an ALPS source checkout or pyalps sdist. "
+      "Configure with -DPYALPS_BUILD_APPLICATIONS=OFF for the core-only package.")
+  endif()
+
+  set(_dmft "${_alps_source_root}/applications/dmft/qmc")
+
+  nanobind_add_module(maxent_c NB_STATIC
+    "${_alps_source_root}/tool/maxent.cpp"
+    "${_alps_source_root}/tool/maxent_helper.cpp"
+    "${_alps_source_root}/tool/maxent_simulation.cpp"
+    "${_alps_source_root}/tool/maxent_parms.cpp")
+
+  nanobind_add_module(cthyb NB_STATIC
+    "${_dmft}/hybridization/hybmain.cpp"
+    "${_dmft}/hybridization/hybsim.cpp"
+    "${_dmft}/hybridization/hyblocal.cpp"
+    "${_dmft}/hybridization/hybint.cpp"
+    "${_dmft}/hybridization/hybfun.cpp"
+    "${_dmft}/hybridization/hybretintfun.cpp"
+    "${_dmft}/hybridization/hybmatrix.cpp"
+    "${_dmft}/hybridization/hybmatrix_ft.cpp"
+    "${_dmft}/hybridization/hybconfig.cpp"
+    "${_dmft}/hybridization/hybupdates.cpp"
+    "${_dmft}/hybridization/hybevaluate.cpp"
+    "${_dmft}/hybridization/hybmeasurements.cpp")
+
+  nanobind_add_module(ctint NB_STATIC
+    "${_dmft}/interaction_expansion2/main.cpp"
+    "${_dmft}/fouriertransform.C"
+    "${_dmft}/interaction_expansion2/auxiliary.cpp"
+    "${_dmft}/interaction_expansion2/observables.cpp"
+    "${_dmft}/interaction_expansion2/fastupdate.cpp"
+    "${_dmft}/interaction_expansion2/selfenergy.cpp"
+    "${_dmft}/interaction_expansion2/solver.cpp"
+    "${_dmft}/interaction_expansion2/io.cpp"
+    "${_dmft}/interaction_expansion2/splines.cpp"
+    "${_dmft}/interaction_expansion2/interaction_expansion.cpp"
+    "${_dmft}/interaction_expansion2/measurements.cpp"
+    "${_dmft}/interaction_expansion2/model.cpp")
+
+  nanobind_add_module(dwa_c NB_STATIC "${_bindings}/apps/dwa.cpp")
+
+  list(APPEND _pyalps_targets maxent_c cthyb ctint dwa_c)
+  foreach(_target IN ITEMS maxent_c cthyb ctint)
+    target_compile_definitions(${_target} PRIVATE BUILD_PYTHON_MODULE)
+    target_include_directories(${_target} PRIVATE "${_bindings}" "${_dmft}")
+  endforeach()
+  target_include_directories(dwa_c PRIVATE "${_alps_source_root}/applications/qmc/dwa")
+endif()
+
+foreach(_target IN LISTS _pyalps_targets)
+  target_include_directories(${_target} PRIVATE
+    ${ALPS_INCLUDE_DIRS}
+    ${ALPS_EXTRA_INCLUDE_DIRS})
+  target_link_libraries(${_target} PRIVATE ${_pyalps_link_libraries})
+  target_link_options(${_target} PRIVATE ${_pyalps_link_options})
+  set_target_properties(${_target} PROPERTIES
+    INSTALL_RPATH "${ALPS_LIBRARY_DIRS};${_alps_hdf5_prefix}/lib")
+endforeach()
+
+set(_extension_dir "pyalps/_ext")
+install(TARGETS ${_pyalps_targets}
+  LIBRARY DESTINATION "${_extension_dir}"
+  RUNTIME DESTINATION "${_extension_dir}")
+
+install(DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}/src/pyalps/"
+  DESTINATION pyalps
+  FILES_MATCHING PATTERN "*.py")
diff --git a/bindings/python/pyalps/README.md b/bindings/python/pyalps/README.md
new file mode 100644
index 000000000..dbde7cdf7
--- /dev/null
+++ b/bindings/python/pyalps/README.md
@@ -0,0 +1,37 @@
+# pyalps
+
+Python applications and libraries for the Algorithms and Libraries for
+Physics Simulations (ALPS) project. Binary wheels are available from PyPI:
+
+```sh
+python -m pip install pyalps
+```
+
+Install `pyalps[plot]` to use the Matplotlib plotting helpers.
+
+The bindings are built as a standalone `scikit-build-core` project using
+nanobind. A source build requires Python 3.10 or newer, CMake 3.18 or newer,
+a C++17 compiler, BLAS/LAPACK, HDF5, and an installed ALPS C++ SDK. Point
+`ALPS_DIR` at the SDK's `share/alps` package directory.
+
+From the repository root:
+
+```sh
+cmake -S . -B _build/alps -G Ninja \
+  -DCMAKE_INSTALL_PREFIX="$PWD/_build/install" \
+  -DCMAKE_CXX_COMPILER_LAUNCHER=ccache \
+  -DALPS_ENABLE_MPI=OFF \
+  -DALPS_BUILD_LIBS_ONLY=ON
+cmake --build _build/alps --target install
+
+ALPS_DIR="$PWD/_build/install/share/alps" \
+  CMAKE_ARGS="-DCMAKE_CXX_COMPILER_LAUNCHER=ccache" \
+  python -m build --wheel bindings/python/pyalps
+```
+
+The wheel is written to `bindings/python/pyalps/dist` and can be installed
+with `python -m pip install`.
+
+`PYALPS_BUILD_APPLICATIONS=ON` is the default and preserves the MaxEnt,
+DWA, CT-HYB, and CT-INT extension modules. Set it to `OFF` through CMake
+configuration for a smaller core-only developer build.
diff --git a/bindings/python/pyalps/cpp/apps/dwa.cpp b/bindings/python/pyalps/cpp/apps/dwa.cpp
new file mode 100644
index 000000000..9016fd6e8
--- /dev/null
+++ b/bindings/python/pyalps/cpp/apps/dwa.cpp
@@ -0,0 +1,109 @@
+/*****************************************************************************
+*
+* ALPS Project Applications: Directed Worm Algorithm
+*
+* Copyright (C) 2013 by Matthias Troyer  <troyer@phys.ethz.ch> ,
+*                       Lode Pollet      <pollet@phys.ethz.ch> ,
+*                       Ping Nang Ma     <pingnang@phys.ethz.ch>
+*               2026 by the ALPS collaboration
+*
+* Permission is hereby granted, free of charge, to any person obtaining
+* a copy of this software and associated documentation files (the "Software"),
+* to deal in the Software without restriction, including without limitation
+* the rights to use, copy, modify, merge, publish, distribute, sublicense,
+* and/or sell copies of the Software, and to permit persons to whom the
+* Software is furnished to do so, subject to the following conditions:
+*
+* The above copyright notice and this permission notice shall be included
+* in all copies or substantial portions of the Software.
+*
+* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+* DEALINGS IN THE SOFTWARE.
+*
+*****************************************************************************/
+// dwa_c — nanobind port.
+//
+// The old boost::python binding registered each std::vector<T> we
+// handed out as a distinct nb::class_<std::vector<T>> carrying the
+// vector_indexing_suite. nanobind's built-in STL caster auto-
+// converts std::vector<T> ↔ Python list for us, so those
+// registrations are retired; any code that wrote
+// `dwa_c.std_vector_double(...)` now just passes / receives a Python
+// list directly.
+#include <nanobind/nanobind.h>
+#include <nanobind/stl/string.h>
+#include <nanobind/stl/vector.h>
+namespace nb = nanobind;
+#include "worldlines.hpp"
+#include "bandstructure.hpp"
+#include <string>
+#include <vector>
+NB_MODULE(dwa_c, m) {
+    m.doc() = "ALPS DWA (directed worm algorithm) Python bindings.";
+    nb::class_<kink>(m, "kink")
+        .def(nb::init<unsigned int>(), nb::arg("siteindicator"))
+        .def(nb::init<unsigned int, double, unsigned short>(),
+             nb::arg("siteindicator"), nb::arg("time"), nb::arg("state"))
+        .def("__repr__",     &kink::representation)
+        .def("siteindicator",&kink::siteindicator)
+        .def("time",         &kink::time)
+        .def("state",        &kink::state);
+    nb::class_<worldlines::location_type>(m, "location_type");
+    nb::class_<worldlines>(m, "worldlines")
+        .def(nb::init<>())
+        .def(nb::init<unsigned int>(), nb::arg("num_sites"))
+        .def("__repr__", &worldlines::representation)
+        .def("load", static_cast<void (worldlines::*)(std::string const &)>(&worldlines::load))
+        .def("save", static_cast<void (worldlines::*)(std::string const &) const>(&worldlines::save))
+        .def("open_worldlines", &worldlines::open_worldlines)
+        .def("worldlines_siteindicator", &worldlines::worldlines_siteindicator)
+        .def("worldlines_time",          &worldlines::worldlines_time)
+        .def("worldlines_state",         &worldlines::worldlines_state)
+        .def("num_sites", &worldlines::num_sites)
+        .def("num_kinks", &worldlines::num_kinks)
+        .def("states",   &worldlines::states)
+        .def("location", &worldlines::location)
+        .def("state_before", &worldlines::state_before)
+        .def("state",        &worldlines::state)
+        .def("is_valid", static_cast<bool (worldlines::*)(unsigned short)>(&worldlines::is_valid));
+    nb::class_<wormpair>(m, "wormpair")
+        .def(nb::init<>())
+        .def(nb::init<worldlines::location_type, kink, bool, bool>())
+        .def("__repr__", &wormpair::representation)
+        .def("wormhead", &wormpair::wormhead)
+        .def("wormtail", &wormpair::wormtail)
+        .def("wormhead_site",    &wormpair::site)
+        .def("wormhead_time",    &wormpair::time)
+        .def("wormhead_forward", &wormpair::forward)
+        .def("wormtail_site", &wormpair::wormtail_site)
+        .def("wormtail_time", &wormpair::wormtail_time)
+        .def("next_partnersite", &wormpair::next_partnersite)
+        .def("next_time",        &wormpair::next_time)
+        .def("wormhead_turns_around",                         &wormpair::wormhead_turns_around)
+        .def("wormhead_moves_to_new_time",                    &wormpair::wormhead_moves_to_new_time)
+        .def("wormhead_inserts_vertex_and_jumps_to_new_site", &wormpair::wormhead_inserts_vertex_and_jumps_to_new_site)
+        .def("wormhead_deletes_vertex_and_jumps_to_new_site", &wormpair::wormhead_deletes_vertex_and_jumps_to_new_site)
+        .def("wormhead_relinks_vertex_and_jumps_to_new_site", &wormpair::wormhead_relinks_vertex_and_jumps_to_new_site)
+        .def("wormhead_crosses_vertex",                       &wormpair::wormhead_crosses_vertex)
+        .def("wormhead_annihilates_wormtail",                 &wormpair::wormhead_annihilates_wormtail);
+    nb::class_<bandstructure>(m, "bandstructure")
+        .def(nb::init<double, double, double, double, unsigned int>(),
+             nb::arg("V0"), nb::arg("lambda"), nb::arg("a"), nb::arg("m"), nb::arg("L"))
+        .def(nb::init<std::vector<double> const &, std::vector<double> const &,
+                      double, double, unsigned int>(),
+             nb::arg("V0"), nb::arg("lambda"), nb::arg("a"), nb::arg("m"), nb::arg("L"))
+        .def("__repr__", static_cast<std::string (bandstructure::*)()>(&bandstructure::representation))
+        .def("t",       static_cast<std::vector<double> (bandstructure::*)()>(&bandstructure::get_t))
+        .def("U",       static_cast<double              (bandstructure::*)()>(&bandstructure::get_U))
+        .def("Ut",      static_cast<std::vector<double> (bandstructure::*)()>(&bandstructure::get_Ut))
+        .def("norm",    static_cast<std::vector<double> (bandstructure::*)()>(&bandstructure::get_norm))
+        .def("q",       static_cast<std::vector<double> (bandstructure::*)(unsigned int)>(&bandstructure::get_q))
+        .def("wk2",     static_cast<std::vector<double> (bandstructure::*)(unsigned int)>(&bandstructure::get_wk2))
+        .def("wk2_c",   static_cast<double              (bandstructure::*)()>(&bandstructure::get_wk2_c))
+        .def("wk2_d",   static_cast<double              (bandstructure::*)()>(&bandstructure::get_wk2_d));
+}
diff --git a/bindings/python/pyalps/cpp/dict_to_params.hpp b/bindings/python/pyalps/cpp/dict_to_params.hpp
new file mode 100644
index 000000000..d50c84e8d
--- /dev/null
+++ b/bindings/python/pyalps/cpp/dict_to_params.hpp
@@ -0,0 +1,35 @@
+// Copyright (C) 2026 by the ALPS collaboration
+// Part of the ALPS Project — see LICENSE.txt for full license text.
+// SPDX-License-Identifier: MIT
+#ifndef PYALPS_DICT_TO_PARAMS_HPP
+#define PYALPS_DICT_TO_PARAMS_HPP
+#include <alps/ngs/params.hpp>
+#include <nanobind/nanobind.h>
+#include <nanobind/stl/string.h>
+#include <nanobind/stl/vector.h>
+#include <string>
+#include <vector>
+namespace pyalps {
+namespace nb = nanobind;
+inline alps::params params_from_dict(nb::dict const & values) {
+    alps::params result;
+    for (auto item : values) {
+        std::string key = nb::cast<std::string>(nb::str(item.first));
+        nb::handle value = item.second;
+        if (nb::isinstance<nb::bool_>(value))
+            result[key] = nb::cast<bool>(value);
+        else if (nb::isinstance<nb::int_>(value))
+            result[key] = nb::cast<long>(value);
+        else if (nb::isinstance<nb::float_>(value))
+            result[key] = nb::cast<double>(value);
+        else if (nb::isinstance<nb::str>(value))
+            result[key] = nb::cast<std::string>(value);
+        else if (nb::isinstance<nb::list>(value) || nb::isinstance<nb::tuple>(value))
+            result[key] = nb::cast<std::vector<double>>(value);
+        else
+            throw nb::type_error(("unsupported parameter type for '" + key + "'").c_str());
+    }
+    return result;
+}
+} // namespace pyalps
+#endif
diff --git a/bindings/python/pyalps/cpp/ngs/accumulator.cpp b/bindings/python/pyalps/cpp/ngs/accumulator.cpp
new file mode 100644
index 000000000..f3ca7f005
--- /dev/null
+++ b/bindings/python/pyalps/cpp/ngs/accumulator.cpp
@@ -0,0 +1,123 @@
+// Copyright (C) 2010 - 2011 by Lukas Gamper <gamperl@gmail.com>
+//               2026 by the ALPS collaboration
+// Part of the ALPS Project — see LICENSE.txt for full license text.
+// SPDX-License-Identifier: MIT
+#include <alps/ngs/accumulator/feature.hpp>
+#include <alps/ngs/accumulator/feature/count.hpp>
+#include <alps/ngs/accumulator/feature/mean.hpp>
+#include <alps/ngs/accumulator/feature/error.hpp>
+#include <alps/ngs/accumulator/feature/binning_analysis.hpp>
+#include <alps/ngs/accumulator/feature/max_num_binning.hpp>
+#include <alps/hdf5/vector.hpp>
+#include <nanobind/nanobind.h>
+#include <sstream>
+#include <string>
+namespace nb = nanobind;
+namespace {
+template <typename T>
+std::string print_value(T const & value) {
+    std::stringstream stream;
+    value.print(stream);
+    return stream.str();
+}
+template <typename Accumulator>
+typename Accumulator::result_type make_result(Accumulator const & value) {
+    return typename Accumulator::result_type(value);
+}
+template <typename Result>
+void bind_result_operators(nb::class_<Result> & cls) {
+    cls
+        .def("__neg__", [](Result value) { value.negate(); return value; })
+        .def("__iadd__", [](Result & self, Result const & other) -> Result & { self += other; return self; }, nb::rv_policy::reference_internal)
+        .def("__iadd__", [](Result & self, double value) -> Result & { self += value; return self; }, nb::rv_policy::reference_internal)
+        .def("__isub__", [](Result & self, Result const & other) -> Result & { self -= other; return self; }, nb::rv_policy::reference_internal)
+        .def("__isub__", [](Result & self, double value) -> Result & { self -= value; return self; }, nb::rv_policy::reference_internal)
+        .def("__imul__", [](Result & self, Result const & other) -> Result & { self *= other; return self; }, nb::rv_policy::reference_internal)
+        .def("__imul__", [](Result & self, double value) -> Result & { self *= value; return self; }, nb::rv_policy::reference_internal)
+        .def("__itruediv__", [](Result & self, Result const & other) -> Result & { self /= other; return self; }, nb::rv_policy::reference_internal)
+        .def("__itruediv__", [](Result & self, double value) -> Result & { self /= value; return self; }, nb::rv_policy::reference_internal)
+        .def("__add__", [](Result value, Result const & other) { value += other; return value; }, nb::is_operator())
+        .def("__add__", [](Result value, double other) { value += other; return value; }, nb::is_operator())
+        .def("__radd__", [](Result value, double other) { value += other; return value; }, nb::is_operator())
+        .def("__sub__", [](Result value, Result const & other) { value -= other; return value; }, nb::is_operator())
+        .def("__sub__", [](Result value, double other) { value -= other; return value; }, nb::is_operator())
+        .def("__rsub__", [](Result value, double other) { value.negate(); value += other; return value; }, nb::is_operator())
+        .def("__mul__", [](Result value, Result const & other) { value *= other; return value; }, nb::is_operator())
+        .def("__mul__", [](Result value, double other) { value *= other; return value; }, nb::is_operator())
+        .def("__rmul__", [](Result value, double other) { value *= other; return value; }, nb::is_operator())
+        .def("__truediv__", [](Result value, Result const & other) { value /= other; return value; }, nb::is_operator())
+        .def("__truediv__", [](Result value, double other) { value /= other; return value; }, nb::is_operator())
+        .def("__rtruediv__", [](Result value, double other) { value.inverse(); value *= other; return value; }, nb::is_operator())
+        .def("sin",  [](Result value) { value.sin();  return value; })
+        .def("cos",  [](Result value) { value.cos();  return value; })
+        .def("tan",  [](Result value) { value.tan();  return value; })
+        .def("sinh", [](Result value) { value.sinh(); return value; })
+        .def("cosh", [](Result value) { value.cosh(); return value; })
+        .def("tanh", [](Result value) { value.tanh(); return value; })
+        .def("abs",  [](Result value) { value.abs();  return value; })
+        .def("sqrt", [](Result value) { value.sqrt(); return value; })
+        .def("log",  [](Result value) { value.log();  return value; });
+}
+template <typename T>
+void bind_serializable(nb::class_<T> & cls) {
+    cls.def("__str__", &print_value<T>)
+       .def("save", &T::save)
+       .def("load", &T::load)
+       .def("reset", &T::reset);
+}
+} // namespace
+NB_MODULE(pyngsaccumulator_c, m) {
+    using namespace alps::accumulator::impl;
+    using count_accumulator = Accumulator<double, alps::accumulator::count_tag, AccumulatorBase<double>>;
+    using count_result = count_accumulator::result_type;
+    nb::class_<count_accumulator> count_acc(m, "count_accumulator");
+    count_acc.def(nb::init<>()).def("__call__", [](count_accumulator & self, double value) { self(value); })
+             .def("result", &make_result<count_accumulator>).def("count", &count_accumulator::count);
+    bind_serializable(count_acc);
+    nb::class_<count_result> count_res(m, "count_result");
+    count_res.def(nb::init<>()).def("count", &count_result::count);
+    bind_serializable(count_res); bind_result_operators(count_res);
+    using mean_accumulator = Accumulator<double, alps::accumulator::mean_tag, count_accumulator>;
+    using mean_result = mean_accumulator::result_type;
+    nb::class_<mean_accumulator> mean_acc(m, "mean_accumulator");
+    mean_acc.def(nb::init<>()).def("__call__", [](mean_accumulator & self, double value) { self(value); })
+            .def("result", &make_result<mean_accumulator>).def("count", &mean_accumulator::count)
+            .def("mean", &mean_accumulator::mean);
+    bind_serializable(mean_acc);
+    nb::class_<mean_result> mean_res(m, "mean_result");
+    mean_res.def(nb::init<>()).def("count", &mean_result::count).def("mean", &mean_result::mean);
+    bind_serializable(mean_res); bind_result_operators(mean_res);
+    using error_accumulator = Accumulator<double, alps::accumulator::error_tag, mean_accumulator>;
+    using error_result = error_accumulator::result_type;
+    nb::class_<error_accumulator> error_acc(m, "error_accumulator");
+    error_acc.def(nb::init<>()).def("__call__", [](error_accumulator & self, double value) { self(value); })
+             .def("result", &make_result<error_accumulator>).def("count", &error_accumulator::count)
+             .def("mean", &error_accumulator::mean).def("error", &error_accumulator::error);
+    bind_serializable(error_acc);
+    nb::class_<error_result> error_res(m, "error_result");
+    error_res.def(nb::init<>()).def("count", &error_result::count).def("mean", &error_result::mean)
+             .def("error", &error_result::error);
+    bind_serializable(error_res); bind_result_operators(error_res);
+    using binning_accumulator = Accumulator<double, alps::accumulator::binning_analysis_tag, error_accumulator>;
+    using binning_result = binning_accumulator::result_type;
+    nb::class_<binning_accumulator> binning_acc(m, "binning_analysis_accumulator");
+    binning_acc.def(nb::init<>()).def("__call__", [](binning_accumulator & self, double value) { self(value); })
+               .def("result", &make_result<binning_accumulator>).def("count", &binning_accumulator::count)
+               .def("mean", &binning_accumulator::mean).def("error", &binning_accumulator::error);
+    bind_serializable(binning_acc);
+    nb::class_<binning_result> binning_res(m, "binning_analysis_result");
+    binning_res.def(nb::init<>()).def("count", &binning_result::count).def("mean", &binning_result::mean)
+               .def("error", &binning_result::error);
+    bind_serializable(binning_res); bind_result_operators(binning_res);
+    using maxbin_accumulator = Accumulator<double, alps::accumulator::max_num_binning_tag, error_accumulator>;
+    using maxbin_result = maxbin_accumulator::result_type;
+    nb::class_<maxbin_accumulator> maxbin_acc(m, "max_num_binning_accumulator");
+    maxbin_acc.def(nb::init<>()).def("__call__", [](maxbin_accumulator & self, double value) { self(value); })
+              .def("result", &make_result<maxbin_accumulator>).def("count", &maxbin_accumulator::count)
+              .def("mean", &maxbin_accumulator::mean).def("error", &maxbin_accumulator::error);
+    bind_serializable(maxbin_acc);
+    nb::class_<maxbin_result> maxbin_res(m, "max_num_binning_result");
+    maxbin_res.def(nb::init<>()).def("count", &maxbin_result::count).def("mean", &maxbin_result::mean)
+              .def("error", &maxbin_result::error);
+    bind_serializable(maxbin_res); bind_result_operators(maxbin_res);
+}
diff --git a/bindings/python/pyalps/cpp/ngs/api.cpp b/bindings/python/pyalps/cpp/ngs/api.cpp
new file mode 100644
index 000000000..09a90dc4b
--- /dev/null
+++ b/bindings/python/pyalps/cpp/ngs/api.cpp
@@ -0,0 +1,32 @@
+// Copyright (C) 2010 - 2011 by Lukas Gamper <gamperl@gmail.com>
+//               Matthias Troyer <troyer@comp-phys.org>
+//               2026       by the ALPS collaboration
+// Part of the ALPS Project — see LICENSE.txt for full license text.
+// SPDX-License-Identifier: MIT
+// <alps/ngs.hpp> umbrella retired in Phase 4 Slice 6 of the NGS
+// retirement (ngs-retirement-scoping.md). This binding's surviving
+// surface (`saveResults`) only needs `alps::mcresults`,
+// `alps::params`, and `alps::hdf5::archive` — pull the narrow
+// headers directly.
+#include <nanobind/nanobind.h>
+#include <nanobind/stl/string.h>
+#include <alps/ngs/mcresults.hpp>
+#include <alps/ngs/params.hpp>
+#include <alps/hdf5/archive.hpp>
+#include <alps/mcbase.hpp>
+namespace nb = nanobind;
+namespace alps {
+    namespace detail {
+        void save_results_export(mcresults const & res, params const & par, alps::hdf5::archive & ar, std::string const & path) {
+            ar["/parameters"] << par;
+            if (res.size())
+                ar[path] << res;
+        }
+    }
+}
+NB_MODULE(pyngsapi_c, m) {
+    m.def("collectResults", [](alps::mcbase const & sim) {
+        return alps::collect_results<alps::mcbase>(sim);
+    });
+    m.def("saveResults", &alps::detail::save_results_export);
+}
diff --git a/bindings/python/pyalps/cpp/ngs/extract_from_pyobject.hpp b/bindings/python/pyalps/cpp/ngs/extract_from_pyobject.hpp
new file mode 100644
index 000000000..3023d4bad
--- /dev/null
+++ b/bindings/python/pyalps/cpp/ngs/extract_from_pyobject.hpp
@@ -0,0 +1,130 @@
+// Copyright (C) 2026 by the ALPS collaboration
+// Part of the ALPS Project — see LICENSE.txt for full license text.
+// SPDX-License-Identifier: MIT
+/// Header-only template: dispatches a nb::object to a visitor based
+/// on the object's Python type name. For numpy arrays the buffer
+/// protocol (PEP 3118) provides the raw data + shape + format string;
+/// for numpy scalars nb::cast handles the conversion through the
+/// scalar's __int__/__float__/__complex__ methods. No dependence on
+/// <numpy/arrayobject.h>.
+#ifndef PYALPS_NGS_EXTRACT_FROM_PYOBJECT_HPP
+#define PYALPS_NGS_EXTRACT_FROM_PYOBJECT_HPP
+    #include <alps/ngs/cast.hpp>
+    #include <alps/ngs/config.hpp>
+    #include <nanobind/nanobind.h>
+    #include <nanobind/ndarray.h>
+    #include <nanobind/stl/complex.h>
+    #include <nanobind/stl/string.h>
+    #include <complex>
+    #include <cstdint>
+    #include <stdexcept>
+    #include <string>
+    #include <vector>
+    namespace alps {
+        namespace detail {
+            namespace nb_ = nanobind;
+            // True iff the ndarray view is C-contiguous: strides match
+            // the canonical row-major layout (rightmost stride =
+            // itemsize, each leftward stride = previous * shape[i+1]).
+            // Zero-rank scalars are trivially contiguous.
+            template <typename Arr>
+            inline bool ndarray_is_c_contiguous(Arr const & arr) {
+                if (arr.ndim() == 0) return true;
+                int64_t expected = 1;
+                for (int64_t i = arr.ndim() - 1; i >= 0; --i) {
+                    if (arr.shape(i) == 0) return true; // empty array
+                    if (arr.stride(i) != expected) return false;
+                    expected *= arr.shape(i);
+                }
+                return true;
+            }
+            /// Dispatches `data` to `visitor` based on Python's type
+            /// name. Visitor must be callable with bool / int / long /
+            /// double / std::complex<double> / std::string /
+            /// nb::list / nb::dict, plus the two-arg numpy form
+            /// `visitor(T const*, std::vector<std::size_t>)` for each
+            /// supported native numpy element type.
+            template<typename T> void extract_from_pyobject_py11(T & visitor, nb_::handle data) {
+                std::string dtype = data.ptr()->ob_type->tp_name;
+                if (dtype == "bool") visitor(nb_::cast<bool>(data));
+                else if (dtype == "int") visitor(nb_::cast<int>(data));
+                else if (dtype == "long") visitor(nb_::cast<long>(data));
+                else if (dtype == "float") visitor(nb_::cast<double>(data));
+                else if (dtype == "complex") visitor(nb_::cast<std::complex<double>>(data));
+                else if (dtype == "str") visitor(nb_::cast<std::string>(data));
+                else if (dtype == "list") visitor(nb_::borrow<nb_::list>(data));
+                else if (dtype == "tuple") {
+                    // materialise the tuple as a list so the visitor only
+                    // needs one sequence overload.
+                    nb_::list as_list = nb_::steal<nb_::list>(
+                        PySequence_List(data.ptr()));
+                    visitor(as_list);
+                }
+                else if (dtype == "dict") visitor(nb_::borrow<nb_::dict>(data));
+                // numpy scalars: extract through the scalar's own
+                // __int__/__float__/__complex__ — no numpy C macros.
+                else if (dtype == "numpy.str_" || dtype == "numpy.str")
+                    visitor(std::string(nb_::cast<std::string>(nb_::str(data.attr("__str__")()))));
+                else if (dtype == "numpy.bool_" || dtype == "numpy.bool")
+                    visitor(nb_::cast<bool>(data));
+                else if (dtype == "numpy.int8")  visitor(static_cast<std::int8_t>(nb_::cast<long>(data)));
+                else if (dtype == "numpy.int16") visitor(static_cast<std::int16_t>(nb_::cast<long>(data)));
+                else if (dtype == "numpy.int32") visitor(static_cast<std::int32_t>(nb_::cast<long>(data)));
+                else if (dtype == "numpy.int64") visitor(static_cast<std::int64_t>(nb_::cast<long long>(data)));
+                else if (dtype == "numpy.uint8")  visitor(static_cast<std::uint8_t>(nb_::cast<unsigned long>(data)));
+                else if (dtype == "numpy.uint16") visitor(static_cast<std::uint16_t>(nb_::cast<unsigned long>(data)));
+                else if (dtype == "numpy.uint32") visitor(static_cast<std::uint32_t>(nb_::cast<unsigned long>(data)));
+                else if (dtype == "numpy.uint64") visitor(static_cast<std::uint64_t>(nb_::cast<unsigned long long>(data)));
+                else if (dtype == "numpy.float32") visitor(static_cast<float>(nb_::cast<double>(data)));
+                else if (dtype == "numpy.float64") visitor(nb_::cast<double>(data));
+                else if (dtype == "numpy.complex64")
+                    visitor(std::complex<float>(
+                          nb_::cast<double>(data.attr("real").attr("__float__")())
+                        , nb_::cast<double>(data.attr("imag").attr("__float__")())
+                    ));
+                else if (dtype == "numpy.complex128")
+                    visitor(nb_::cast<std::complex<double>>(data));
+                else if (dtype == "numpy.ndarray") {
+                    // Raw buffer access via nb::ndarray, with a strict
+                    // dtype match — nb::cast<nb::ndarray<T>>(arr) of a
+                    // mismatched-dtype array silently coerces (e.g.
+                    // int → bool yields all-true), so we inspect
+                    // .dtype() ourselves and pick the matching arm.
+                    // We require C-contiguity; the typical save path
+                    // is bulk contiguous data and silently copying
+                    // behind the user's back was the old
+                    // PyArray_GETCONTIGUOUS behaviour we don't want
+                    // to inherit.
+                    auto arr_any = nb_::cast<nb_::ndarray<nb_::ro, nb_::c_contig>>(data);
+                    std::vector<std::size_t> sizes;
+                    sizes.reserve(arr_any.ndim());
+                    for (std::size_t i = 0; i < arr_any.ndim(); ++i)
+                        sizes.push_back(static_cast<std::size_t>(arr_any.shape(i)));
+                    auto dt = arr_any.dtype();
+                    #define DISPATCH_DTYPE(T)                                                           \
+                        if (dt == nb_::dtype<T>())                                                      \
+                            return visitor(static_cast<T const *>(arr_any.data()), sizes);
+                    DISPATCH_DTYPE(bool)
+                    DISPATCH_DTYPE(signed char)
+                    DISPATCH_DTYPE(unsigned char)
+                    DISPATCH_DTYPE(short)
+                    DISPATCH_DTYPE(unsigned short)
+                    DISPATCH_DTYPE(int)
+                    DISPATCH_DTYPE(unsigned)
+                    DISPATCH_DTYPE(long)
+                    DISPATCH_DTYPE(unsigned long)
+                    DISPATCH_DTYPE(long long)
+                    DISPATCH_DTYPE(unsigned long long)
+                    DISPATCH_DTYPE(float)
+                    DISPATCH_DTYPE(double)
+                    DISPATCH_DTYPE(std::complex<float>)
+                    DISPATCH_DTYPE(std::complex<double>)
+                    #undef DISPATCH_DTYPE
+                    throw std::runtime_error(
+                        "Unknown numpy element dtype at save site" + ALPS_STACKTRACE);
+                } else
+                    throw std::runtime_error("Unsupported type: " + dtype + ALPS_STACKTRACE);
+            }
+        } // namespace detail
+    } // namespace alps
+#endif  // PYALPS_NGS_EXTRACT_FROM_PYOBJECT_HPP
diff --git a/bindings/python/pyalps/cpp/ngs/hdf5.cpp b/bindings/python/pyalps/cpp/ngs/hdf5.cpp
new file mode 100644
index 000000000..05eeefe92
--- /dev/null
+++ b/bindings/python/pyalps/cpp/ngs/hdf5.cpp
@@ -0,0 +1,317 @@
+// Copyright (C) 2010 - 2012 by Lukas Gamper <gamperl@gmail.com>
+//               Matthias Troyer <troyer@comp-phys.org>
+//               2026       by the ALPS collaboration
+// Part of the ALPS Project — see LICENSE.txt for full license text.
+// SPDX-License-Identifier: MIT
+// Save path: extract_from_pyobject_py11 dispatches a nb::handle to a
+// visitor that writes a concrete C++ value. Load path: dispatches on
+// the archive's inspected type and reads into a concrete C++ type
+// before wrapping it back as a nb::object.
+//
+// Exception translation: pyalps/hdf5.py creates ArchiveError etc. and
+// calls register_archive_exception_type(id, type); the translators
+// below fire PyErr_SetString against whichever Python type was handed
+// in.
+#include <nanobind/nanobind.h>
+#include <nanobind/stl/string.h>
+#include <nanobind/stl/vector.h>
+#include <alps/hdf5/archive.hpp>
+#include <alps/hdf5/pair.hpp>
+#include <alps/hdf5/pointer.hpp>
+#include <alps/hdf5/vector.hpp>
+#include <alps/hdf5/complex.hpp>
+#include <alps/ngs/stacktrace.hpp>
+#include "extract_from_pyobject.hpp"
+#include "../numpy_compat.hpp"
+#include <array>
+#include <complex>
+#include <cstddef>
+#include <cstdint>
+#include <stdexcept>
+#include <string>
+#include <vector>
+namespace nb = nanobind;
+namespace alps {
+    namespace detail {
+        // Save-side visitor: receives a concrete C++ value (or a
+        // nb::list / nb::dict) from extract_from_pyobject_py11 and
+        // writes it to the archive at `path`.
+        struct hdf5_save_py11_visitor {
+            alps::hdf5::archive & ar;
+            std::string path;
+            template <typename U>
+            void operator()(U const & v) const {
+                ar[path] << v;
+            }
+            template <typename U>
+            void operator()(U const * ptr, std::vector<std::size_t> const & sizes) const {
+                // Use make_pvp(path, ptr, size-vector) to preserve the
+                // dimensional shape — a plain vector<U> flatten would
+                // round-trip the data but lose the rank.
+                ar << alps::make_pvp(path, ptr, sizes);
+            }
+            void operator()(nb::list const & l) const {
+                // Order: flat numeric first, then nested numeric, then
+                // strings. Heterogeneous / deeply-nested / mixed-type
+                // lists fall through to the descent branch below which
+                // stores each entry under a numeric child path.
+                try { ar[path] << nb::cast<std::vector<double>>(l); return; }
+                catch (nb::cast_error const &) {}
+                try { ar[path] << nb::cast<std::vector<int>>(l); return; }
+                catch (nb::cast_error const &) {}
+                try { ar[path] << nb::cast<std::vector<std::complex<double>>>(l); return; }
+                catch (nb::cast_error const &) {}
+                try { ar[path] << nb::cast<std::vector<std::vector<double>>>(l); return; }
+                catch (nb::cast_error const &) {}
+                try { ar[path] << nb::cast<std::vector<std::vector<int>>>(l); return; }
+                catch (nb::cast_error const &) {}
+                try { ar[path] << nb::cast<std::vector<std::string>>(l); return; }
+                catch (nb::cast_error const &) {}
+                // Inhomogeneous — recurse per-element into
+                // <path>/<index>, letting each entry be stored as its
+                // own native type.
+                ar.create_group(path);
+                Py_ssize_t i = 0;
+                for (auto item : l) {
+                    std::string child = path + "/" + std::to_string(static_cast<long long>(i++));
+                    hdf5_save_py11_visitor child_visitor{ar, child};
+                    extract_from_pyobject_py11(child_visitor, item);
+                }
+            }
+            void operator()(nb::dict const & d) const {
+                // Store a dict as a group with one child per key. Keys
+                // are stringified (HDF5 paths are strings), values go
+                // through the same save dispatch recursively.
+                ar.create_group(path);
+                for (auto item : d) {
+                    std::string key = nb::cast<std::string>(nb::str(item.first));
+                    std::string child = path + "/" + key;
+                    hdf5_save_py11_visitor child_visitor{ar, child};
+                    extract_from_pyobject_py11(child_visitor, item.second);
+                }
+            }
+        };
+        std::string python_hdf5_get_filename(alps::hdf5::archive & ar) {
+            return ar.get_filename();
+        }
+        void python_hdf5_save(alps::hdf5::archive & ar,
+                              std::string const & path,
+                              nb::handle data) {
+            hdf5_save_py11_visitor visitor{ar, path};
+            extract_from_pyobject_py11(visitor, data);
+        }
+        // Helper: load a multi-dim HDF5 dataset of element type T
+        // into a flat std::vector, then wrap as a numpy array with
+        // the original shape so that Python sees a 2-D np.array for
+        // rank-2 writes etc. Preserves the dimensionality encoded on
+        // the save path (alps::make_pvp(path, ptr, size-vector)).
+        template <typename T>
+        nb::object load_nd_array(alps::hdf5::archive & ar,
+                                 std::string const & path,
+                                 std::vector<std::size_t> const & shape) {
+            std::size_t total = 1;
+            for (auto s : shape) total *= s;
+            std::vector<T> flat(total);
+            if (shape.size() <= 1) {
+                // vector<T> overload works directly.
+                ar[path] >> flat;
+            } else {
+                // make_pvp with explicit size-vector to read a
+                // multi-dim dataset into a flat buffer.
+                ar >> alps::make_pvp(path, flat.data(), shape);
+            }
+            return alps::python::make_numpy_array<T>(flat.data(), shape);
+        }
+        nb::object python_hdf5_load_impl(alps::hdf5::archive & ar,
+                                         std::string const & path);
+        nb::object python_hdf5_load(alps::hdf5::archive & ar,
+                                    std::string const & path) {
+            return python_hdf5_load_impl(ar, path);
+        }
+        nb::object python_hdf5_load_impl(alps::hdf5::archive & ar,
+                                         std::string const & path) {
+            // Groups (not datasets) get loaded recursively. Children
+            // whose names are consecutive decimal integers starting at 0
+            // are recovered as a Python list (preserving round-trip for
+            // list-saved-as-group); otherwise a dict.
+            if (ar.is_group(path)) {
+                auto children = ar.list_children(path);
+                bool list_shaped = true;
+                for (std::size_t i = 0; list_shaped && i < children.size(); ++i) {
+                    if (children[i] != std::to_string(i))
+                        list_shaped = false;
+                }
+                if (list_shaped) {
+                    nb::list result;
+                    for (auto const & child : children)
+                        result.append(
+                            python_hdf5_load_impl(ar, path + "/" + child));
+                    return nb::object(std::move(result));
+                } else {
+                    nb::dict result;
+                    for (auto const & child : children)
+                        result[nb::str(child.c_str())] =
+                            python_hdf5_load_impl(ar, path + "/" + child);
+                    return nb::object(std::move(result));
+                }
+            }
+            // Complex values have a quirky HDF5 representation: a
+            // single complex is stored as rank-1 dims=[2] (real,imag)
+            // and a 2x2 array of complex as rank-3 dims=[2,2,2]. So
+            // is_scalar returns false for a scalar complex — branch
+            // on is_complex first and use the rank minus 1 (stripping
+            // the trailing complex-pair dim) to tell scalar from
+            // array.
+            if (ar.is_complex(path)) {
+                auto ext = ar.extent(path);
+                if (ext.size() == 1) {
+                    std::complex<double> v; ar[path] >> v; return nb::cast(v);
+                }
+                std::vector<std::size_t> shape(ext.begin(), ext.end() - 1);
+                return load_nd_array<std::complex<double>>(ar, path, shape);
+            }
+            // Convenience macros for the scalar path: check each
+            // candidate integer width in turn (numpy's default int is
+            // platform-dependent — int64 on macOS/Linux, int32 on
+            // Windows — so we can't rely on just `int` matching).
+            #define TRY_SCALAR(T)                                                                \
+                if (ar.is_datatype<T>(path)) { T v; ar[path] >> v; return nb::cast(v); }
+            if (ar.is_scalar(path)) {
+                TRY_SCALAR(std::string)
+                TRY_SCALAR(double)
+                TRY_SCALAR(float)
+                TRY_SCALAR(bool)
+                TRY_SCALAR(std::int64_t)
+                TRY_SCALAR(std::int32_t)
+                TRY_SCALAR(std::int16_t)
+                TRY_SCALAR(std::int8_t)
+                TRY_SCALAR(std::uint64_t)
+                TRY_SCALAR(std::uint32_t)
+                TRY_SCALAR(std::uint16_t)
+                TRY_SCALAR(std::uint8_t)
+                throw std::runtime_error(
+                    "Unknown HDF5 scalar type at " + path + ALPS_STACKTRACE);
+            } else {
+                // String datasets don't map to nb::ndarray<T> the way
+                // numeric types do; keep the flat vector path.
+                if (ar.is_datatype<std::string>(path)) {
+                    std::vector<std::string> v; ar[path] >> v; return nb::cast(v);
+                }
+                auto shape = ar.extent(path);
+                if (ar.is_datatype<double>(path))        return load_nd_array<double>(ar, path, shape);
+                if (ar.is_datatype<float>(path))         return load_nd_array<float>(ar, path, shape);
+                if (ar.is_datatype<std::int64_t>(path))  return load_nd_array<std::int64_t>(ar, path, shape);
+                if (ar.is_datatype<std::int32_t>(path))  return load_nd_array<std::int32_t>(ar, path, shape);
+                if (ar.is_datatype<std::int16_t>(path))  return load_nd_array<std::int16_t>(ar, path, shape);
+                if (ar.is_datatype<std::int8_t>(path))   return load_nd_array<std::int8_t>(ar, path, shape);
+                if (ar.is_datatype<std::uint64_t>(path)) return load_nd_array<std::uint64_t>(ar, path, shape);
+                if (ar.is_datatype<std::uint32_t>(path)) return load_nd_array<std::uint32_t>(ar, path, shape);
+                if (ar.is_datatype<std::uint16_t>(path)) return load_nd_array<std::uint16_t>(ar, path, shape);
+                if (ar.is_datatype<std::uint8_t>(path))  return load_nd_array<std::uint8_t>(ar, path, shape);
+                throw std::runtime_error(
+                    "Unknown HDF5 vector type at " + path + ALPS_STACKTRACE);
+            }
+            #undef TRY_SCALAR
+        }
+        nb::list python_hdf5_extent(alps::hdf5::archive & ar,
+                                    std::string const & path) {
+            nb::list result;
+            std::vector<std::size_t> ext = ar.extent(path);
+            if (ar.is_complex(path)) {
+                if (ext.size() > 1)
+                    ext.pop_back();
+                else
+                    ext.back() = 1;
+            }
+            for (auto const & s : ext)
+                result.append(s);
+            return result;
+        }
+        // Python exception types registered by pyalps.hdf5 at import
+        // time. Translators below fire PyErr_SetString against these
+        // pre-registered PyObject*'s so the Python side sees its own
+        // subclasses (ArchiveError / ArchiveNotFound / ...).
+        std::array<PyObject *, 6> exception_type = {};
+        #define TRANSLATE_CPP_ERROR_TO_PYTHON(T, ID)                            \
+            static void translate_ ## T (hdf5:: T const & e) {                  \
+                std::string message =                                           \
+                    std::string(e.what()).substr(                               \
+                        0, std::string(e.what()).find_first_of('\n'));          \
+                PyErr_SetString(exception_type[ID] ? exception_type[ID]         \
+                                                   : PyExc_RuntimeError,       \
+                                message.c_str());                               \
+            }
+        TRANSLATE_CPP_ERROR_TO_PYTHON(archive_error, 0)
+        TRANSLATE_CPP_ERROR_TO_PYTHON(archive_not_found, 1)
+        TRANSLATE_CPP_ERROR_TO_PYTHON(archive_closed, 2)
+        TRANSLATE_CPP_ERROR_TO_PYTHON(invalid_path, 3)
+        TRANSLATE_CPP_ERROR_TO_PYTHON(path_not_found, 4)
+        TRANSLATE_CPP_ERROR_TO_PYTHON(wrong_type, 5)
+        #undef TRANSLATE_CPP_ERROR_TO_PYTHON
+        void register_exception_type(int id, nb::object type) {
+            if (id < 0 || id >= static_cast<int>(exception_type.size()))
+                throw std::out_of_range(
+                    "register_archive_exception_type: id out of range");
+            // Py_INCREF the incoming type so it survives past this call
+            // (we're keeping a raw PyObject* in a static array).
+            Py_INCREF(type.ptr());
+            exception_type[id] = type.ptr();
+        }
+    }
+}
+NB_MODULE(pyngshdf5_c, m) {
+    // Install the six C++→Python exception translators. Each calls the
+    // matching translate_* above, which forwards to whichever Python
+    // class was registered via register_archive_exception_type. If
+    // pyalps/hdf5.py hasn't run yet, the translator falls back to
+    // RuntimeError so the module is safely loadable on its own.
+    nb::register_exception_translator(
+        [](const std::exception_ptr &p, void * /*payload*/) {
+            try { std::rethrow_exception(p); }
+            catch (alps::hdf5::archive_not_found const & e) {
+                alps::detail::translate_archive_not_found(e);
+            } catch (alps::hdf5::archive_closed const & e) {
+                alps::detail::translate_archive_closed(e);
+            } catch (alps::hdf5::invalid_path const & e) {
+                alps::detail::translate_invalid_path(e);
+            } catch (alps::hdf5::path_not_found const & e) {
+                alps::detail::translate_path_not_found(e);
+            } catch (alps::hdf5::wrong_type const & e) {
+                alps::detail::translate_wrong_type(e);
+            } catch (alps::hdf5::archive_error const & e) {
+                // Base class — must be caught LAST since the specialized
+                // types above inherit from it.
+                alps::detail::translate_archive_error(e);
+            }
+        });
+    m.def("register_archive_exception_type",
+          &alps::detail::register_exception_type);
+    nb::class_<alps::hdf5::archive>(m, "hdf5_archive_impl")
+        .def(nb::init<std::string, std::string>())
+        .def("__deepcopy__",
+             // copy.deepcopy() hands us (self, memo); memo unused.
+             [](alps::hdf5::archive const & self, nb::handle /*memo*/) {
+                 return alps::hdf5::archive(self);
+             })
+        .def_prop_ro("filename", &alps::detail::python_hdf5_get_filename)
+        .def_prop_ro("context",  &alps::hdf5::archive::get_context)
+        .def_prop_ro("is_open",  &alps::hdf5::archive::is_open)
+        .def("set_context",     &alps::hdf5::archive::set_context)
+        .def("is_group",        &alps::hdf5::archive::is_group)
+        .def("is_data",         &alps::hdf5::archive::is_data)
+        .def("is_attribute",    &alps::hdf5::archive::is_attribute)
+        .def("close",           &alps::hdf5::archive::close)
+        .def("extent",          &alps::detail::python_hdf5_extent)
+        .def("dimensions",      &alps::hdf5::archive::dimensions)
+        .def("is_scalar",       &alps::hdf5::archive::is_scalar)
+        .def("is_complex",      &alps::hdf5::archive::is_complex)
+        .def("is_null",         &alps::hdf5::archive::is_null)
+        .def("list_children",   &alps::hdf5::archive::list_children)
+        .def("list_attributes", &alps::hdf5::archive::list_attributes)
+        .def("__setitem__",     &alps::detail::python_hdf5_save)
+        .def("__getitem__",     &alps::detail::python_hdf5_load)
+        .def("create_group",    &alps::hdf5::archive::create_group)
+        .def("delete_data",     &alps::hdf5::archive::delete_data)
+        .def("delete_group",    &alps::hdf5::archive::delete_group)
+        .def("delete_attribute",&alps::hdf5::archive::delete_attribute);
+}
diff --git a/bindings/python/pyalps/cpp/ngs/mcbase.cpp b/bindings/python/pyalps/cpp/ngs/mcbase.cpp
new file mode 100644
index 000000000..0c4b704c7
--- /dev/null
+++ b/bindings/python/pyalps/cpp/ngs/mcbase.cpp
@@ -0,0 +1,173 @@
+/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+ *                                                                                 *
+ * ALPS Project: Algorithms and Libraries for Physics Simulations                  *
+ *                                                                                 *
+ * ALPS Libraries                                                                  *
+ *                                                                                 *
+ * Copyright (C) 2010 - 2011 by Lukas Gamper <gamperl@gmail.com>                   *
+ *                              Matthias Troyer <troyer@comp-phys.org>             *
+ *               2026       by the ALPS collaboration                              *
+ *                                                                                 *
+ * Permission is hereby granted, free of charge, to any person obtaining           *
+ * a copy of this software and associated documentation files (the "Software"),    *
+ * to deal in the Software without restriction, including without limitation       *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,        *
+ * and/or sell copies of the Software, and to permit persons to whom the           *
+ * Software is furnished to do so, subject to the following conditions:            *
+ *                                                                                 *
+ * The above copyright notice and this permission notice shall be included         *
+ * in all copies or substantial portions of the Software.                          *
+ *                                                                                 *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS         *
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,     *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE     *
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER          *
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING         *
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER             *
+ * DEALINGS IN THE SOFTWARE.                                                       *
+ *                                                                                 *
+ * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
+// pyngsbase_c — nanobind port.
+//
+// Trampoline (PyMCBase) forwards the three pure-virtual mcbase methods
+// (update / measure / fraction_completed) back into the Python subclass
+// through nanobind's trampoline support. The old wrapper<mcbase>
+// pattern becomes a standard trampoline-plus-alias pair.
+//
+// Params ingestion: the public alps::mcbase ctor wants an alps::params.
+// libalps still declares a params(boost::python::dict) ctor in its
+// header, but we don't want to drag boost::python through the
+// nanobind bindings. Instead, we convert nb::dict → alps::params at the
+// binding boundary by iterating and setitem-ing concrete C++ values
+// (int/float/bool/str/list). That sidesteps the cross-registry issue
+// and keeps the libalps ABI untouched.
+#define PY_ARRAY_UNIQUE_SYMBOL pyngsbase_PyArrayHandle
+#include <alps/mcbase.hpp>
+#include <alps/hdf5/archive.hpp>
+#include <nanobind/nanobind.h>
+#include <nanobind/make_iterator.h>
+#include <nanobind/stl/function.h>
+#include <nanobind/stl/string.h>
+#include <nanobind/stl/vector.h>
+#include <nanobind/trampoline.h>
+namespace nb = nanobind;
+#ifdef ALPS_HAVE_MPI
+    #include <boost/mpi.hpp>
+#endif
+#include <cstddef>
+#include <functional>
+#include <stdexcept>
+#include <string>
+#include <vector>
+namespace alps {
+    namespace detail {
+        // Convert a Python dict into an alps::params, extracting concrete
+        // C++ values for each entry. This mirrors what the libalps
+        // params(boost::python::dict) ctor does, but without routing the
+        // nb::object through the boost::python::object variant alternative
+        // — everything stays within the nanobind type registry.
+        inline alps::params py_dict_to_params(nb::dict const & d) {
+            alps::params p;
+            for (auto item : d) {
+                std::string k = nb::cast<std::string>(nb::str(item.first));
+                nb::handle v = item.second;
+                if (nb::isinstance<nb::bool_>(v))
+                    p[k] = nb::cast<bool>(v);
+                else if (nb::isinstance<nb::int_>(v))
+                    p[k] = nb::cast<int>(v);
+                else if (nb::isinstance<nb::float_>(v))
+                    p[k] = nb::cast<double>(v);
+                else if (nb::isinstance<nb::str>(v))
+                    p[k] = nb::cast<std::string>(v);
+                else if (nb::isinstance<nb::list>(v) || nb::isinstance<nb::tuple>(v))
+                    p[k] = nb::cast<std::vector<double>>(v);
+                else
+                    throw nb::type_error((
+                        "unsupported type for key '" + k +
+                        "' in params dict (expected bool/int/float/str/list)").c_str());
+            }
+            return p;
+        }
+    }
+}
+namespace alps {
+    // Trampoline: holds Python overrides for pure-virtuals. The
+    // protected mcbase members (random / parameters / measurements)
+    // are accessed via lambdas in the binding below, which friend-in
+    // through PyMCBase (a protected member is visible to a derived
+    // class's own member functions / friends).
+    class PyMCBase : public mcbase {
+        public:
+            NB_TRAMPOLINE(mcbase, 3);
+            #ifdef ALPS_HAVE_MPI
+                PyMCBase(nb::dict const & arg,
+                         std::size_t seed_offset = 42,
+                         boost::mpi::communicator const & /*comm*/ = boost::mpi::communicator())
+                    : mcbase(alps::detail::py_dict_to_params(arg), seed_offset)
+                {}
+            #else
+                PyMCBase(nb::dict const & arg, std::size_t seed_offset = 42)
+                    : mcbase(alps::detail::py_dict_to_params(arg), seed_offset)
+                {}
+            #endif
+            void update() override {
+                NB_OVERRIDE_PURE(update);
+            }
+            void measure() override {
+                NB_OVERRIDE_PURE(measure);
+            }
+            double fraction_completed() const override {
+                NB_OVERRIDE_PURE(fraction_completed);
+            }
+            // Accessors for protected mcbase members. Called from the
+            // binding lambdas below (they friend-in through PyMCBase).
+            alps::random01 & get_random() { return random; }
+            mcbase::parameters_type & get_parameters() { return parameters; }
+            alps::mcobservables & get_measurements() { return measurements; }
+            // mcbase::run takes a std::function<bool()>; wrap a Python
+            // callable so the stop_callback can be driven from Python.
+            bool run_py(nb::object stop_callback) {
+                return mcbase::run([stop_callback]() -> bool {
+                    nb::gil_scoped_acquire gil;
+                    return nb::cast<bool>(stop_callback());
+                });
+            }
+    };
+}
+NB_MODULE(pyngsbase_c, m) {
+    nb::class_<alps::mcbase>(m, "_mcbase", nb::never_destruct());
+    nb::class_<alps::PyMCBase, alps::mcbase>(m, "mcbase")
+        // Always expose the (dict, seed_offset) form from Python. When
+        // ALPS_HAVE_MPI is on we'd *like* to offer an optional
+        // communicator too, but boost::mpi::communicator is not a
+        // nanobind-registered type so nb::arg(..).default_value() can't
+        // materialise it. MPI simulations that actually need to hand
+        // Python a communicator should do so from C++ using the
+        // extended trampoline ctor directly.
+        .def(nb::init<nb::dict const &, std::size_t>(),
+             nb::arg("dict"),
+             nb::arg("seed_offset") = 42)
+        .def_prop_ro(
+            "random",
+            [](alps::PyMCBase & self) -> alps::random01 & { return self.get_random(); },
+            nb::rv_policy::reference_internal)
+        .def_prop_ro(
+            "parameters",
+            [](alps::PyMCBase & self) -> alps::mcbase::parameters_type & { return self.get_parameters(); },
+            nb::rv_policy::reference_internal)
+        .def_prop_ro(
+            "measurements",
+            [](alps::PyMCBase & self) -> alps::mcobservables & { return self.get_measurements(); },
+            nb::rv_policy::reference_internal)
+        .def("run",
+             [](alps::PyMCBase & self, nb::object cb) { return self.run_py(std::move(cb)); })
+        // Pure-virtual methods: bound on the base class; the trampoline's
+        // The trampoline forwards the call into the Python subclass.
+        .def("update",             &alps::mcbase::update)
+        .def("measure",            &alps::mcbase::measure)
+        .def("fraction_completed", &alps::mcbase::fraction_completed)
+        .def("save", static_cast<void (alps::mcbase::*)(alps::hdf5::archive &) const>(
+                         &alps::mcbase::save))
+        .def("load", static_cast<void (alps::mcbase::*)(alps::hdf5::archive &)>(
+                         &alps::mcbase::load));
+}
diff --git a/bindings/python/pyalps/cpp/ngs/observable.cpp b/bindings/python/pyalps/cpp/ngs/observable.cpp
new file mode 100644
index 000000000..dde673f21
--- /dev/null
+++ b/bindings/python/pyalps/cpp/ngs/observable.cpp
@@ -0,0 +1,81 @@
+/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+ *                                                                                 *
+ * ALPS Project: Algorithms and Libraries for Physics Simulations                  *
+ *                                                                                 *
+ * ALPS Libraries                                                                  *
+ *                                                                                 *
+ * Copyright (C) 2010 - 2011 by Lukas Gamper <gamperl@gmail.com>                   *
+ *                              Matthias Troyer <troyer@comp-phys.org>             *
+ *               2026       by the ALPS collaboration                              *
+ *                                                                                 *
+ * Permission is hereby granted, free of charge, to any person obtaining           *
+ * a copy of this software and associated documentation files (the "Software"),    *
+ * to deal in the Software without restriction, including without limitation       *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,        *
+ * and/or sell copies of the Software, and to permit persons to whom the           *
+ * Software is furnished to do so, subject to the following conditions:            *
+ *                                                                                 *
+ * The above copyright notice and this permission notice shall be included         *
+ * in all copies or substantial portions of the Software.                          *
+ *                                                                                 *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS         *
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,     *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE     *
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER          *
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING         *
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER             *
+ * DEALINGS IN THE SOFTWARE.                                                       *
+ *                                                                                 *
+ * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
+// pyngsobservable_c — nanobind port.
+#define PY_ARRAY_UNIQUE_SYMBOL pyngsobservable_PyArrayHandle
+#include <alps/hdf5/archive.hpp>
+#include <alps/hdf5/complex.hpp>
+#include <alps/ngs/mcobservable.hpp>
+#include <nanobind/nanobind.h>
+#include <nanobind/ndarray.h>
+#include <nanobind/stl/string.h>
+namespace nb = nanobind;
+#include <alps/alea/detailedbinning.h>
+#include <memory>
+#include <stdexcept>
+#include <string>
+#include <valarray>
+namespace alps {
+    namespace detail {
+        void observable_append(alps::mcobservable & self, nb::object const & data) {
+            if (nb::isinstance<nb::int_>(data) || nb::isinstance<nb::float_>(data)) {
+                self << nb::cast<double>(data);
+                return;
+            }
+            try {
+                auto values = nb::cast<nb::ndarray<double, nb::c_contig>>(data);
+                self << std::valarray<double>(values.data(), values.size());
+            } catch (nb::cast_error const &) {
+                throw nb::type_error("observable samples must be numeric scalars or contiguous float64 arrays");
+            }
+        }
+        void observable_load(alps::mcobservable & self, alps::hdf5::archive & ar, std::string const & path) {
+            std::string current = ar.get_context();
+            ar.set_context(path);
+            self.load(ar);
+            ar.set_context(current);
+        }
+        alps::mcobservable create_RealObservable_export(std::string name) {
+            return alps::mcobservable(std::make_shared<alps::RealObservable>(name).get());
+        }
+        alps::mcobservable create_RealVectorObservable_export(std::string name) {
+            return alps::mcobservable(std::make_shared<alps::RealVectorObservable>(name).get());
+        }
+    }
+}
+NB_MODULE(pyngsobservable_c, m) {
+    m.def("createRealObservable",       &alps::detail::create_RealObservable_export);
+    m.def("createRealVectorObservable", &alps::detail::create_RealVectorObservable_export);
+    nb::class_<alps::mcobservable>(m, "observable")
+        .def("append",          &alps::detail::observable_append)
+        .def("merge",           &alps::mcobservable::merge)
+        .def("save",            &alps::mcobservable::save)
+        .def("load",            &alps::detail::observable_load)
+        .def("addToObservable", &alps::detail::observable_load);
+}
diff --git a/bindings/python/pyalps/cpp/ngs/observables.cpp b/bindings/python/pyalps/cpp/ngs/observables.cpp
new file mode 100644
index 000000000..7b121fe64
--- /dev/null
+++ b/bindings/python/pyalps/cpp/ngs/observables.cpp
@@ -0,0 +1,108 @@
+/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+ *                                                                                 *
+ * ALPS Project: Algorithms and Libraries for Physics Simulations                  *
+ *                                                                                 *
+ * ALPS Libraries                                                                  *
+ *                                                                                 *
+ * Copyright (C) 2010 - 2011 by Lukas Gamper <gamperl@gmail.com>                   *
+ *                              Matthias Troyer <troyer@comp-phys.org>             *
+ *               2026       by the ALPS collaboration                              *
+ *                                                                                 *
+ * Permission is hereby granted, free of charge, to any person obtaining           *
+ * a copy of this software and associated documentation files (the "Software"),    *
+ * to deal in the Software without restriction, including without limitation       *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,        *
+ * and/or sell copies of the Software, and to permit persons to whom the           *
+ * Software is furnished to do so, subject to the following conditions:            *
+ *                                                                                 *
+ * The above copyright notice and this permission notice shall be included         *
+ * in all copies or substantial portions of the Software.                          *
+ *                                                                                 *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS         *
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,     *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE     *
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER          *
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING         *
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER             *
+ * DEALINGS IN THE SOFTWARE.                                                       *
+ *                                                                                 *
+ * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
+// pyngsobservables_c — nanobind port.
+//
+// Unlike mcresults, alps::mcobservables doesn't override erase() so we
+// can stand on a straight nb::class_<> with hand-written map methods
+// that match the map_indexing_suite surface area. We keep it explicit
+// (rather than nb::bind_map) because the class also carries non-map
+// methods (reset/save/load/__lshift__/create*) that need to live on
+// the same binding, and mixing bind_map with extra defs is noisy.
+#define PY_ARRAY_UNIQUE_SYMBOL pyngsobservables_PyArrayHandle
+#include <alps/hdf5/archive.hpp>
+#include <alps/ngs/mcobservable.hpp>
+#include <alps/ngs/mcobservables.hpp>
+#include <alps/ngs/observablewrappers.hpp>
+#include <nanobind/nanobind.h>
+#include <nanobind/make_iterator.h>
+#include <nanobind/stl/string.h>
+namespace nb = nanobind;
+#include <cstdint>
+#include <string>
+namespace {
+void mcobservables_load(alps::mcobservables & self, alps::hdf5::archive & ar, std::string const & path) {
+    std::string current = ar.get_context();
+    ar.set_context(path);
+    self.load(ar);
+    ar.set_context(current);
+}
+void createRealObservable(alps::mcobservables & self, std::string const & name, std::uint32_t binnum) {
+    self << alps::ngs::RealObservable(name, binnum);
+}
+void createRealVectorObservable(alps::mcobservables & self, std::string const & name, std::uint32_t binnum) {
+    self << alps::ngs::RealVectorObservable(name, binnum);
+}
+void addObservable(alps::mcobservables & self, nb::object const & obj) {
+    // Mirror boost::python::call_method<void>(obj, "addToObservables", ref(self)):
+    // bounce the call back into Python, passing `self` by reference.
+    obj.attr("addToObservables")(nb::cast(&self, nb::rv_policy::reference));
+}
+}  // namespace
+NB_MODULE(pyngsobservables_c, m) {
+    nb::class_<alps::mcobservables>(m, "observables")
+        .def(nb::init<>())
+        .def("__len__",      [](alps::mcobservables const & self) { return self.size(); })
+        .def("__contains__", [](alps::mcobservables const & self, std::string const & k) {
+                                 return self.has(k);
+                             })
+        .def("__getitem__",  [](alps::mcobservables & self, std::string const & k) -> alps::mcobservable & {
+                                 if (!self.has(k))
+                                     throw nb::key_error(k.c_str());
+                                 return self[k];
+                             },
+                             nb::rv_policy::reference_internal)
+        .def("__setitem__",  [](alps::mcobservables & self, std::string const & k, alps::mcobservable const & v) {
+                                 self.insert(k, v);
+                             })
+        .def("__iter__",     [](alps::mcobservables & self) {
+                                 return nb::make_key_iterator(nb::type<alps::mcobservables>(), "key_iterator", self.begin(), self.end());
+                             },
+                             nb::keep_alive<0, 1>())
+        .def("keys",         [](alps::mcobservables & self) {
+                                 return nb::make_key_iterator(nb::type<alps::mcobservables>(), "key_iterator", self.begin(), self.end());
+                             },
+                             nb::keep_alive<0, 1>())
+        .def("values",       [](alps::mcobservables & self) {
+                                 return nb::make_value_iterator(nb::type<alps::mcobservables>(), "value_iterator", self.begin(), self.end());
+                             },
+                             nb::keep_alive<0, 1>())
+        .def("items",        [](alps::mcobservables & self) {
+                                 return nb::make_iterator(nb::type<alps::mcobservables>(), "item_iterator", self.begin(), self.end());
+                             },
+                             nb::keep_alive<0, 1>())
+        .def("reset", &alps::mcobservables::reset, nb::arg("equilibrated") = false)
+        .def("save",  &alps::mcobservables::save)
+        .def("load",  &mcobservables_load)
+        .def("__lshift__", &addObservable)
+        .def("createRealObservable",       &createRealObservable,
+             nb::arg("name"), nb::arg("binnum") = 0)
+        .def("createRealVectorObservable", &createRealVectorObservable,
+             nb::arg("name"), nb::arg("binnum") = 0);
+}
diff --git a/bindings/python/pyalps/cpp/ngs/params.cpp b/bindings/python/pyalps/cpp/ngs/params.cpp
new file mode 100644
index 000000000..f6d9c8957
--- /dev/null
+++ b/bindings/python/pyalps/cpp/ngs/params.cpp
@@ -0,0 +1,160 @@
+// Copyright (C) 2010 - 2011 by Lukas Gamper <gamperl@gmail.com>
+//               Matthias Troyer <troyer@comp-phys.org>
+//               2026       by the ALPS collaboration
+// Part of the ALPS Project — see LICENSE.txt for full license text.
+// SPDX-License-Identifier: MIT
+#include <nanobind/nanobind.h>
+#include <nanobind/make_iterator.h>
+#include <nanobind/stl/string.h>
+#include <nanobind/stl/vector.h>
+#include <alps/hdf5/archive.hpp>
+#include <alps/ngs/params.hpp>
+#include <alps/ngs/detail/paramvalue.hpp>
+#include <complex>
+#include <sstream>
+#include <stdexcept>
+#include <string>
+#include <boost/variant/apply_visitor.hpp>
+#include <boost/variant/static_visitor.hpp>
+#include <vector>
+namespace nb = nanobind;
+namespace {
+// Convert a Python dict into an alps::params. Same shape as the
+// helper in mcbase.cpp but kept local to params.cpp so a change to
+// the dispatch (e.g. adding complex support) can stay in one place
+// alongside the other setitem logic.
+alps::params py_dict_to_params(nb::dict const & d);
+// Walk the paramvalue variant and wrap each native alternative as a
+// nb::object. Called from __getitem__.
+struct paramvalue_to_py_visitor : boost::static_visitor<nb::object> {
+    template <typename T>
+    nb::object operator()(T const & value) const {
+        return nb::cast(value);
+    }
+};
+nb::object paramvalue_to_py(alps::detail::paramvalue const & pv) {
+    return boost::apply_visitor(
+        paramvalue_to_py_visitor(),
+        static_cast<alps::detail::paramvalue_base const &>(pv));
+}
+// Deposit a native C++ value from a Python object into the paramvalue
+// via paramproxy's templated operator=.
+void params_setitem(alps::params & self, nb::object const & key_obj, nb::object const & value) {
+    std::string key = nb::cast<std::string>(nb::str(key_obj));
+    if (nb::isinstance<nb::bool_>(value))
+        self[key] = nb::cast<bool>(value);
+    else if (nb::isinstance<nb::int_>(value))
+        self[key] = nb::cast<int>(value);
+    else if (nb::isinstance<nb::float_>(value))
+        self[key] = nb::cast<double>(value);
+    else if (nb::isinstance<nb::str>(value))
+        self[key] = nb::cast<std::string>(value);
+    else if (nb::isinstance<nb::list>(value) || nb::isinstance<nb::tuple>(value)) {
+        // Heuristic: try doubles first, strings as fallback.
+        try {
+            self[key] = nb::cast<std::vector<double>>(value);
+        } catch (nb::cast_error &) {
+            self[key] = nb::cast<std::vector<std::string>>(value);
+        }
+    } else {
+        throw nb::type_error("unsupported value type for params[]");
+    }
+}
+nb::object params_getitem(alps::params & self, nb::object const & key_obj) {
+    std::string key = nb::cast<std::string>(nb::str(key_obj));
+    if (!self.defined(key))
+        return nb::none();
+    // params doesn't expose the underlying map directly, but
+    // paramiterator yields (key, paramvalue) pairs; walk it to find the
+    // entry and hand the variant to paramvalue_to_py.
+    for (auto it = self.begin(); it != self.end(); ++it)
+        if (it->first == key)
+            return paramvalue_to_py(it->second);
+    return nb::none();  // defensive — defined()==true should guarantee a hit
+}
+void params_delitem(alps::params & self, nb::object const & key_obj) {
+    self.erase(nb::cast<std::string>(nb::str(key_obj)));
+}
+bool params_contains(alps::params & self, nb::object const & key_obj) {
+    return self.defined(nb::cast<std::string>(nb::str(key_obj)));
+}
+nb::object value_or_default(alps::params & self, nb::object const & key, nb::object const & dflt) {
+    return params_contains(self, key) ? params_getitem(self, key) : dflt;
+}
+void params_load(alps::params & self, alps::hdf5::archive & ar, std::string const & path) {
+    std::string current = ar.get_context();
+    ar.set_context(path);
+    self.load(ar);
+    ar.set_context(current);
+}
+std::string params_print(alps::params & self) {
+    std::stringstream ss;
+    ss << self;
+    return ss.str();
+}
+// deepcopy support — nanobind passes (self, memo); memo unused.
+alps::params params_deepcopy(alps::params const & self, nb::handle /*memo*/) {
+    return alps::params(self);
+}
+// Materialise an alps::params from a Python dict. Re-uses the same
+// type dispatch as params_setitem so a round-tripped dict-built
+// params contains exactly the same variant alternatives.
+alps::params py_dict_to_params(nb::dict const & d) {
+    alps::params p;
+    for (auto item : d) {
+        std::string k = nb::cast<std::string>(nb::str(item.first));
+        nb::handle v = item.second;
+        if (nb::isinstance<nb::bool_>(v))
+            p[k] = nb::cast<bool>(v);
+        else if (nb::isinstance<nb::int_>(v))
+            p[k] = nb::cast<int>(v);
+        else if (nb::isinstance<nb::float_>(v))
+            p[k] = nb::cast<double>(v);
+        else if (nb::isinstance<nb::str>(v))
+            p[k] = nb::cast<std::string>(v);
+        else if (nb::isinstance<nb::list>(v) || nb::isinstance<nb::tuple>(v)) {
+            try { p[k] = nb::cast<std::vector<double>>(v); }
+            catch (nb::cast_error &) {
+                p[k] = nb::cast<std::vector<std::string>>(v);
+            }
+        } else {
+            throw nb::type_error(
+                ("unsupported value type for params key '" + k + "'").c_str());
+        }
+    }
+    return p;
+}
+}  // namespace
+NB_MODULE(pyngsparams_c, m) {
+    nb::class_<alps::params>(m, "params")
+        .def(nb::init<>())
+        .def("__init__",
+             [](alps::params * self, nb::dict const & d) {
+                 new (self) alps::params(py_dict_to_params(d));
+             },
+             nb::arg("dict"))
+        .def(nb::init<alps::hdf5::archive, std::string const &>(),
+             nb::arg("archive"),
+             nb::arg("path") = std::string("/parameters"))
+        .def("__len__",      [](alps::params const & self) { return self.size(); })
+        .def("__deepcopy__", &params_deepcopy)
+        .def("__getitem__",  &params_getitem)
+        .def("__setitem__",  &params_setitem)
+        .def("__delitem__",  &params_delitem)
+        .def("__contains__", &params_contains)
+        .def("__iter__",     [](alps::params & self) {
+                                 // paramiterator yields pair<string const, paramvalue>;
+                                 // make_key_iterator projects out pair.first.
+                                 return nb::make_key_iterator(
+                                     nb::type<alps::params>(),
+                                     "key_iterator",
+                                     self.begin(), self.end());
+                             },
+                             nb::keep_alive<0, 1>())
+        .def("__str__",      &params_print)
+        .def("valueOrDefault", &value_or_default)
+        .def("save",         &alps::params::save)
+        .def("load",         &params_load,
+             nb::arg("archive"),
+             nb::arg("path") = std::string("/parameters"));
+}
diff --git a/bindings/python/pyalps/cpp/ngs/random01.cpp b/bindings/python/pyalps/cpp/ngs/random01.cpp
new file mode 100644
index 000000000..4efc47dee
--- /dev/null
+++ b/bindings/python/pyalps/cpp/ngs/random01.cpp
@@ -0,0 +1,24 @@
+// Copyright (C) 2010 - 2013 by Lukas Gamper <gamperl@gmail.com>
+//               Matthias Troyer <troyer@comp-phys.org>
+//               2026       by the ALPS collaboration
+// Part of the ALPS Project — see LICENSE.txt for full license text.
+// SPDX-License-Identifier: MIT
+#include <nanobind/nanobind.h>
+#include <nanobind/stl/string.h>
+#include <alps/ngs/random01.hpp>
+#include <alps/hdf5/archive.hpp>
+namespace nb = nanobind;
+NB_MODULE(pyngsrandom01_c, m) {
+    nb::class_<alps::random01>(m, "random01")
+        .def(nb::init<int>(), nb::arg("seed") = 42)
+        .def("__deepcopy__",
+             // copy.deepcopy() passes (self, memo); memo is unused.
+             [](alps::random01 const & self, nb::handle /*memo*/) {
+                 return alps::random01(self);
+             })
+        .def("__call__",
+             static_cast<alps::random01::result_type (alps::random01::*)()>(
+                 &alps::random01::operator()))
+        .def("save", &alps::random01::save)
+        .def("load", &alps::random01::load);
+}
diff --git a/bindings/python/pyalps/cpp/ngs/result.cpp b/bindings/python/pyalps/cpp/ngs/result.cpp
new file mode 100644
index 000000000..5d7d6248e
--- /dev/null
+++ b/bindings/python/pyalps/cpp/ngs/result.cpp
@@ -0,0 +1,181 @@
+// Copyright (C) 2010 - 2011 by Lukas Gamper <gamperl@gmail.com>
+//               Matthias Troyer <troyer@comp-phys.org>
+//               2026       by the ALPS collaboration
+// Part of the ALPS Project — see LICENSE.txt for full license text.
+// SPDX-License-Identifier: MIT
+#include <nanobind/nanobind.h>
+#include <nanobind/ndarray.h>
+#include <nanobind/stl/string.h>
+#include <alps/hdf5/archive.hpp>
+#include <alps/ngs/cast.hpp>
+#include <alps/ngs/mcresult.hpp>
+#include <alps/ngs/mcobservable.hpp>
+#include <alps/ngs/stacktrace.hpp>
+#include <cmath>
+#include <cstddef>
+#include <cstring>
+#include <sstream>
+#include <stdexcept>
+#include <string>
+#include <vector>
+namespace nb = nanobind;
+namespace {
+// Allocate a heap-owned 1-D numpy array of T and copy `n` elements from `src`.
+template <typename T>
+nb::ndarray<nb::numpy, T> make_1d(T const * src, std::size_t n) {
+    T * buf = new T[n ? n : 1];
+    if (n) std::memcpy(buf, src, n * sizeof(T));
+    nb::capsule owner(buf, [](void * p) noexcept {
+        delete[] static_cast<T *>(p);
+    });
+    std::size_t shape[1] = { n };
+    return nb::ndarray<nb::numpy, T>(buf, 1, shape, owner);
+}
+} // namespace
+namespace alps {
+    namespace detail {
+        template <typename T> std::string short_print_python(T const & value) {
+            return cast<std::string>(value);
+        }
+        template <typename T> std::string short_print_python(std::vector<T> const & value) {
+            switch (value.size()) {
+                case 0:
+                    return "[]";
+                case 1:
+                    return "[" + short_print_python(value.front()) + "]";
+                case 2:
+                    return "[" + short_print_python(value.front()) + "," + short_print_python(value.back()) + "]";
+                default:
+                    return "[" + short_print_python(value.front()) + ",.." + short_print_python(value.size()) + "..," + short_print_python(value.back()) + "]";
+            }
+        }
+        inline nb::object vec_to_numpy(std::vector<double> const & v) {
+            return nb::cast(make_1d<double>(v.data(), v.size()));
+        }
+        std::string mcresult_print(alps::mcresult const & self) {
+            if (self.count() == 0)
+                return "No Measurements";
+            else if (self.is_type<double>())
+                return short_print_python(self.mean<double>()) + "(" + short_print_python(self.count()) + ") "
+                    + "+/-" + short_print_python(self.error<double>()) + " "
+                    + short_print_python(self.bins<double>()) + "#" + short_print_python(self.bin_size());
+            else if (self.is_type<std::vector<double> >())
+                return short_print_python(self.mean<std::vector<double> >()) + "(" + short_print_python(self.count()) + ") "
+                    + "+/-" + short_print_python(self.error<std::vector<double> >()) + " "
+                    + short_print_python(self.bins<std::vector<double> >()) + "#" + short_print_python(self.bin_size());
+            else
+                throw std::runtime_error("Unsupported type." + ALPS_STACKTRACE);
+        }
+        nb::object mcresult_mean(alps::mcresult const & self) {
+            if (self.is_type<double>())
+                return nb::float_(self.mean<double>());
+            else if (self.is_type<std::vector<double> >())
+                return vec_to_numpy(self.mean<std::vector<double> >());
+            else
+                throw std::runtime_error("Unsupported type." + ALPS_STACKTRACE);
+        }
+        nb::object mcresult_error(alps::mcresult const & self) {
+            if (self.is_type<double>())
+                return nb::float_(self.error<double>());
+            else if (self.is_type<std::vector<double> >())
+                return vec_to_numpy(self.error<std::vector<double> >());
+            else
+                throw std::runtime_error("Unsupported type." + ALPS_STACKTRACE);
+        }
+        nb::object mcresult_tau(alps::mcresult const & self) {
+            if (self.is_type<double>())
+                return nb::float_(self.tau<double>());
+            else if (self.is_type<std::vector<double> >())
+                return vec_to_numpy(self.tau<std::vector<double> >());
+            else
+                throw std::runtime_error("Unsupported type." + ALPS_STACKTRACE);
+        }
+        nb::object mcresult_variance(alps::mcresult const & self) {
+            if (self.is_type<double>())
+                return nb::float_(self.variance<double>());
+            else if (self.is_type<std::vector<double> >())
+                return vec_to_numpy(self.variance<std::vector<double> >());
+            else
+                throw std::runtime_error("Unsupported type." + ALPS_STACKTRACE);
+        }
+        nb::object mcresult_bins(alps::mcresult const & self) {
+            if (self.is_type<double>())
+                return vec_to_numpy(self.bins<double>());
+            else
+                throw std::runtime_error("Unsupported type." + ALPS_STACKTRACE);
+        }
+        alps::mcresult observable2result_export(alps::mcobservable const & obs) {
+            return alps::mcresult(obs);
+        }
+    }
+}
+NB_MODULE(pyngsresult_c, m) {
+    using namespace alps;
+    using R = alps::mcresult;
+    m.def("observable2result", &alps::detail::observable2result_export);
+    nb::class_<R>(m, "result")
+        .def(nb::init<>())
+        .def(nb::init<R const &>())
+        .def("__repr__", &alps::detail::mcresult_print)
+        .def("__deepcopy__",
+             [](R const & self, nb::handle /*memo*/) {
+                 return R(self);
+             })
+        .def("__abs__", static_cast<R(*)(R)>(&abs))
+        .def("__pow__", static_cast<R(*)(R, double)>(&pow))
+        .def_prop_ro("mean",     &alps::detail::mcresult_mean)
+        .def_prop_ro("error",    &alps::detail::mcresult_error)
+        .def_prop_ro("tau",      &alps::detail::mcresult_tau)
+        .def_prop_ro("variance", &alps::detail::mcresult_variance)
+        .def_prop_ro("bins",     &alps::detail::mcresult_bins)
+        .def_prop_ro("count",    &R::count)
+        // mcresult's unary +/- operate on non-const self and return a
+        // reference (not a new value). Wrap them in lambdas that
+        // return a fresh copy, which is what Python's +obj/-obj expect.
+        .def("__pos__", [](R self) { return +self; })
+        .def("__neg__", [](R self) { return -self; })
+        // In-place operators — return self by reference so the original
+        // object is modified in place (Python's __i*__ semantics).
+        .def("__iadd__", [](R & self, R const & o) -> R & { return self += o; }, nb::is_operator())
+        .def("__iadd__", [](R & self, double o)    -> R & { return self += o; }, nb::is_operator())
+        .def("__isub__", [](R & self, R const & o) -> R & { return self -= o; }, nb::is_operator())
+        .def("__isub__", [](R & self, double o)    -> R & { return self -= o; }, nb::is_operator())
+        .def("__imul__", [](R & self, R const & o) -> R & { return self *= o; }, nb::is_operator())
+        .def("__imul__", [](R & self, double o)    -> R & { return self *= o; }, nb::is_operator())
+        .def("__itruediv__", [](R & self, R const & o) -> R & { return self /= o; }, nb::is_operator())
+        .def("__itruediv__", [](R & self, double o)    -> R & { return self /= o; }, nb::is_operator())
+        // Binary operators — left and right forms. nb::is_operator()
+        // marks them as Python operator overloads so mixed-type
+        // failures return NotImplemented rather than raising TypeError
+        // (giving Python's reflected operator machinery a chance).
+        .def("__add__",  [](R const & a, R const & b) { return a + b; }, nb::is_operator())
+        .def("__radd__", [](R const & a, R const & b) { return b + a; }, nb::is_operator())
+        .def("__add__",  [](R const & a, double b)    { return a + b; }, nb::is_operator())
+        .def("__radd__", [](R const & a, double b)    { return b + a; }, nb::is_operator())
+        .def("__sub__",  [](R const & a, R const & b) { return a - b; }, nb::is_operator())
+        .def("__rsub__", [](R const & a, R const & b) { return b - a; }, nb::is_operator())
+        .def("__sub__",  [](R const & a, double b)    { return a - b; }, nb::is_operator())
+        .def("__rsub__", [](R const & a, double b)    { return b - a; }, nb::is_operator())
+        .def("__mul__",  [](R const & a, R const & b) { return a * b; }, nb::is_operator())
+        .def("__rmul__", [](R const & a, R const & b) { return b * a; }, nb::is_operator())
+        .def("__mul__",  [](R const & a, double b)    { return a * b; }, nb::is_operator())
+        .def("__rmul__", [](R const & a, double b)    { return b * a; }, nb::is_operator())
+        .def("__truediv__",  [](R const & a, R const & b) { return a / b; }, nb::is_operator())
+        .def("__rtruediv__", [](R const & a, R const & b) { return b / a; }, nb::is_operator())
+        .def("__truediv__",  [](R const & a, double b)    { return a / b; }, nb::is_operator())
+        .def("__rtruediv__", [](R const & a, double b)    { return b / a; }, nb::is_operator())
+        .def("sq",   static_cast<R(*)(R)>(&sq))
+        .def("cb",   static_cast<R(*)(R)>(&cb))
+        .def("sqrt", static_cast<R(*)(R)>(&sqrt))
+        .def("cbrt", static_cast<R(*)(R)>(&cbrt))
+        .def("exp",  static_cast<R(*)(R)>(&exp))
+        .def("log",  static_cast<R(*)(R)>(&log))
+        .def("sin",  static_cast<R(*)(R)>(&sin))
+        .def("cos",  static_cast<R(*)(R)>(&cos))
+        .def("tan",  static_cast<R(*)(R)>(&tan))
+        .def("sinh", static_cast<R(*)(R)>(&sinh))
+        .def("cosh", static_cast<R(*)(R)>(&cosh))
+        .def("tanh", static_cast<R(*)(R)>(&tanh))
+        .def("save", &R::save)
+        .def("load", &R::load);
+}
diff --git a/bindings/python/pyalps/cpp/ngs/results.cpp b/bindings/python/pyalps/cpp/ngs/results.cpp
new file mode 100644
index 000000000..1b4b7d8ad
--- /dev/null
+++ b/bindings/python/pyalps/cpp/ngs/results.cpp
@@ -0,0 +1,86 @@
+// Copyright (C) 2010 - 2011 by Lukas Gamper <gamperl@gmail.com>
+//               Matthias Troyer <troyer@comp-phys.org>
+//               2026       by the ALPS collaboration
+// Part of the ALPS Project — see LICENSE.txt for full license text.
+// SPDX-License-Identifier: MIT
+// A generic map binder cannot be used for alps::mcresults because
+// mcresults::erase(std::string const &) shadows the std::map::erase(iterator)
+// that bind_map relies on for __delitem__. Synthesise the dict-like surface
+// by hand instead. (Same applies to nanobind's bind_map.)
+#define PY_ARRAY_UNIQUE_SYMBOL pyngsresults_PyArrayHandle
+#include <nanobind/nanobind.h>
+#include <nanobind/make_iterator.h>
+#include <nanobind/stl/string.h>
+#include <alps/hdf5.hpp>
+#include <alps/ngs/mcresults.hpp>
+#include <sstream>
+#include <stdexcept>
+#include <string>
+namespace nb = nanobind;
+namespace alps {
+    namespace detail {
+        std::string mcresults_print(alps::mcresults & self) {
+            std::stringstream sstr;
+            sstr << self;
+            return sstr.str();
+        }
+        void mcresults_load(alps::mcresults & self, alps::hdf5::archive & ar, std::string const & path) {
+            std::string current = ar.get_context();
+            ar.set_context(path);
+            self.load(ar);
+            ar.set_context(current);
+        }
+    }
+}
+NB_MODULE(pyngsresults_c, m) {
+    nb::class_<alps::mcresults>(m, "results")
+        .def("__len__",      [](alps::mcresults const & self) { return self.size(); })
+        .def("__contains__", [](alps::mcresults const & self, std::string const & k) {
+                                 return self.has(k);
+                             })
+        .def("__getitem__",  [](alps::mcresults & self, std::string const & k) -> alps::mcresult const & {
+                                 if (!self.has(k))
+                                     throw nb::key_error(k.c_str());
+                                 return self[k];
+                             },
+                             nb::rv_policy::reference_internal)
+        .def("__setitem__",  [](alps::mcresults & self, std::string const & k, alps::mcresult const & v) {
+                                 self.insert(k, v);
+                             })
+        .def("__delitem__",  [](alps::mcresults & self, std::string const & k) {
+                                 if (!self.has(k))
+                                     throw nb::key_error(k.c_str());
+                                 self.erase(k);
+                             })
+        .def("__iter__",     [](alps::mcresults & self) {
+                                 return nb::make_key_iterator(
+                                     nb::type<alps::mcresults>(),
+                                     "key_iterator",
+                                     self.begin(), self.end());
+                             },
+                             nb::keep_alive<0, 1>())
+        .def("keys",         [](alps::mcresults & self) {
+                                 return nb::make_key_iterator(
+                                     nb::type<alps::mcresults>(),
+                                     "key_iterator",
+                                     self.begin(), self.end());
+                             },
+                             nb::keep_alive<0, 1>())
+        .def("values",       [](alps::mcresults & self) {
+                                 return nb::make_value_iterator(
+                                     nb::type<alps::mcresults>(),
+                                     "value_iterator",
+                                     self.begin(), self.end());
+                             },
+                             nb::keep_alive<0, 1>())
+        .def("items",        [](alps::mcresults & self) {
+                                 return nb::make_iterator(
+                                     nb::type<alps::mcresults>(),
+                                     "item_iterator",
+                                     self.begin(), self.end());
+                             },
+                             nb::keep_alive<0, 1>())
+        .def("__str__",      &alps::detail::mcresults_print)
+        .def("save",         &alps::mcresults::save)
+        .def("load",         &alps::detail::mcresults_load);
+}
diff --git a/bindings/python/pyalps/cpp/numpy_compat.hpp b/bindings/python/pyalps/cpp/numpy_compat.hpp
new file mode 100644
index 000000000..b23cafc0e
--- /dev/null
+++ b/bindings/python/pyalps/cpp/numpy_compat.hpp
@@ -0,0 +1,95 @@
+// Copyright (C) 2026 by the ALPS collaboration
+// Part of the ALPS Project — see LICENSE.txt for full license text.
+// SPDX-License-Identifier: MIT
+//
+// Numpy interop without numpy headers. Construct numpy.ndarray
+// instances from C++ buffers and consume incoming numpy arrays through
+// nb::ndarray's DLPack/buffer view. The numpy package itself is loaded
+// at runtime via nb::module_::import_("numpy"); pyalps already requires
+// numpy as a runtime dependency.
+#ifndef ALPS_PYTHON_NUMPY_COMPAT_HPP
+#define ALPS_PYTHON_NUMPY_COMPAT_HPP
+#include <nanobind/nanobind.h>
+#include <nanobind/ndarray.h>
+#include <nanobind/stl/string.h>
+#include <complex>
+#include <cstddef>
+#include <cstdint>
+#include <cstring>
+#include <utility>
+#include <vector>
+namespace alps {
+    namespace python {
+        namespace nb_ = nanobind;
+        // numpy dtype strings, indexed by the corresponding C++ type.
+        // Used by make_numpy_array() / as_contiguous() to drive the
+        // numpy.empty(dtype=…) / numpy.ascontiguousarray(dtype=…) calls.
+        template <typename T> struct numpy_dtype;
+        template <> struct numpy_dtype<bool>                 { static constexpr char const* name = "bool"; };
+        template <> struct numpy_dtype<std::int8_t>          { static constexpr char const* name = "int8"; };
+        template <> struct numpy_dtype<std::int16_t>         { static constexpr char const* name = "int16"; };
+        template <> struct numpy_dtype<std::int32_t>         { static constexpr char const* name = "int32"; };
+        template <> struct numpy_dtype<std::int64_t>         { static constexpr char const* name = "int64"; };
+        template <> struct numpy_dtype<std::uint8_t>         { static constexpr char const* name = "uint8"; };
+        template <> struct numpy_dtype<std::uint16_t>        { static constexpr char const* name = "uint16"; };
+        template <> struct numpy_dtype<std::uint32_t>        { static constexpr char const* name = "uint32"; };
+        template <> struct numpy_dtype<std::uint64_t>        { static constexpr char const* name = "uint64"; };
+        template <> struct numpy_dtype<float>                { static constexpr char const* name = "float32"; };
+        template <> struct numpy_dtype<double>               { static constexpr char const* name = "float64"; };
+        template <> struct numpy_dtype<std::complex<float>>  { static constexpr char const* name = "complex64"; };
+        template <> struct numpy_dtype<std::complex<double>> { static constexpr char const* name = "complex128"; };
+        // Allocates numpy.empty(shape, dtype=numpy_dtype<T>::name) and
+        // memcpy's `data` (length = product(shape)) into it. Returns
+        // a writable numpy.ndarray.
+        template <typename T>
+        inline nb_::object make_numpy_array(T const* data,
+                                            std::vector<std::size_t> const& shape) {
+            nb_::object np = nb_::module_::import_("numpy");
+            nb_::tuple shape_tuple = nb_::steal<nb_::tuple>(PyTuple_New(static_cast<Py_ssize_t>(shape.size())));
+            for (std::size_t i = 0; i < shape.size(); ++i)
+                PyTuple_SET_ITEM(shape_tuple.ptr(), static_cast<Py_ssize_t>(i),
+                                 PyLong_FromUnsignedLongLong(shape[i]));
+            nb_::object arr = np.attr("empty")(
+                shape_tuple, nb_::arg("dtype") = numpy_dtype<T>::name);
+            // Bridge the freshly-allocated numpy buffer through nb::ndarray
+            // to get a writable raw pointer.
+            auto nd = nb_::cast<nb_::ndarray<T, nb_::c_contig>>(arr);
+            std::size_t total = 1;
+            for (auto s : shape) total *= s;
+            if (total > 0)
+                std::memcpy(nd.data(), data, total * sizeof(T));
+            return arr;
+        }
+        template <typename T>
+        inline nb_::object make_numpy_array(std::vector<T> const& v) {
+            return make_numpy_array<T>(v.data(), {v.size()});
+        }
+        // Strong-ref'd C-contiguous view onto a numpy array of dtype T.
+        // The owner handle keeps the array alive for the lifetime of
+        // the view; data() / shape() / ndim() forward to the ndarray.
+        template <typename T>
+        struct contiguous_view {
+            nb_::object owner;
+            nb_::ndarray<T, nb_::c_contig> nd;
+            T const* data() const { return nd.data(); }
+            std::size_t ndim() const { return nd.ndim(); }
+            std::size_t shape(int i) const { return nd.shape(i); }
+        };
+        // Coerces `obj` to a C-contiguous numpy.ndarray of dtype T via
+        // numpy.ascontiguousarray. Always produces a contiguous +
+        // correctly-typed buffer (numpy copies if the input doesn't
+        // already match). Equivalent in spirit to nanobind's
+        // py::array_t<T, py::array::c_style | py::array::forcecast>
+        // parameter form, just routed through numpy at runtime instead
+        // of through the numpy C headers at compile time.
+        template <typename T>
+        inline contiguous_view<T> as_contiguous(nb_::handle obj) {
+            nb_::object np = nb_::module_::import_("numpy");
+            nb_::object arr = np.attr("ascontiguousarray")(
+                obj, nb_::arg("dtype") = numpy_dtype<T>::name);
+            auto nd = nb_::cast<nb_::ndarray<T, nb_::c_contig>>(arr);
+            return contiguous_view<T>{std::move(arr), std::move(nd)};
+        }
+    } // namespace python
+} // namespace alps
+#endif // ALPS_PYTHON_NUMPY_COMPAT_HPP
diff --git a/bindings/python/pyalps/cpp/pyalea.cpp b/bindings/python/pyalps/cpp/pyalea.cpp
new file mode 100644
index 000000000..b5ecd5607
--- /dev/null
+++ b/bindings/python/pyalps/cpp/pyalea.cpp
@@ -0,0 +1,368 @@
+// Copyright (C) 1994-2010 by Ping Nang Ma <pingnang@itp.phys.ethz.ch>,
+//               Lukas Gamper <gamperl@gmail.com>,
+//               Matthias Troyer <troyer@itp.phys.ethz.ch>,
+//               Maximilian Poprawe <poprawem@ethz.ch>
+//               2026       by the ALPS collaboration
+// Part of the ALPS Project — see LICENSE.txt for full license text.
+// SPDX-License-Identifier: MIT
+#include <nanobind/nanobind.h>
+#include <nanobind/ndarray.h>
+#include <nanobind/stl/pair.h>
+#include <nanobind/stl/string.h>
+#include <nanobind/stl/vector.h>
+#include <alps/alea/detailedbinning.h>
+#include <alps/alea/mcanalyze.hpp>
+#include <alps/alea/mcdata.hpp>
+#include <alps/alea/value_with_error.hpp>
+#include <alps/hdf5.hpp>
+#include <alps/numeric/vector_functions.hpp>
+#include "numpy_compat.hpp"
+#include "save_observable_to_hdf5.hpp"
+#include <alps/random.h>
+#include <cstddef>
+#include <sstream>
+#include <stdexcept>
+#include <string>
+#include <utility>
+#include <valarray>
+#include <vector>
+namespace nb = nanobind;
+namespace alps {
+    namespace alea {
+        // Wraps a scalar-valarray observable (RealVectorObservable /
+        // RealVectorTimeSeriesObservable) so Python sees numpy arrays
+        // in / out instead of std::valarray<double>.
+        template<typename T>
+        class WrappedValarrayObservable {
+            using element_type = typename T::value_type::value_type;
+          public:
+            WrappedValarrayObservable(std::string const & name, int s = 0)
+                : obs(name, s)
+            {}
+            // Copy the ndarray into a valarray and feed it to the observable.
+            void push(nb::handle arr) {
+                auto view = alps::python::as_contiguous<double>(arr);
+                if (view.ndim() != 1)
+                    throw std::invalid_argument(
+                        "RealVectorObservable.push: expected 1-D array");
+                std::size_t n = static_cast<std::size_t>(view.shape(0));
+                std::valarray<element_type> v(n);
+                double const * data = view.data();
+                for (std::size_t i = 0; i < n; ++i)
+                    v[i] = static_cast<element_type>(data[i]);
+                obs << v;
+            }
+            std::string representation() const { return obs.representation(); }
+            // Turn an alps::numeric std::valarray-like view into a
+            // 1-D numpy.ndarray (dtype=float64) by copying through
+            // numpy.empty + buffer protocol.
+            template <typename U>
+            static nb::object _to_numpy(U const & v) {
+                std::size_t n = static_cast<std::size_t>(v.size());
+                std::vector<double> tmp(n);
+                for (std::size_t i = 0; i < n; ++i)
+                    tmp[i] = static_cast<double>(v[i]);
+                return alps::python::make_numpy_array<double>(tmp.data(), {n});
+            }
+            nb::object mean()     const { return _to_numpy(obs.mean());     }
+            nb::object error()    const { return _to_numpy(obs.error());    }
+            nb::object tau()      const { return _to_numpy(obs.tau());      }
+            nb::object variance() const { return _to_numpy(obs.variance()); }
+            void save(std::string const & filename) const {
+                alps::hdf5::archive ar(filename, "a");
+                ar["/simulation/results/" + obs.representation()] << obs;
+            }
+            typename T::count_type       count()            const { return obs.count(); }
+            typename T::convergence_type converged_errors() const { return obs.converged_errors(); }
+          private:
+            T obs;
+        };
+    } // namespace alea
+} // namespace alps
+namespace {
+// Build a 1-D numpy.ndarray (dtype=float64) from any sequence-like
+// alps container (std::vector<double>, std::valarray<double>).
+template <typename Container>
+nb::object seq_to_numpy(Container const & v) {
+    std::size_t n = static_cast<std::size_t>(v.size());
+    std::vector<double> tmp(n);
+    for (std::size_t i = 0; i < n; ++i)
+        tmp[i] = static_cast<double>(v[i]);
+    return alps::python::make_numpy_array<double>(tmp.data(), {n});
+}
+// Copy a numpy array into a std::vector<ValueType>. Used when the
+// caller still instantiates mctimeseries<ValueType> from a
+// Python array.
+template <typename ValueType>
+std::vector<ValueType>
+numpy_to_vector(nb::handle arr) {
+    auto view = alps::python::as_contiguous<double>(arr);
+    if (view.ndim() != 1)
+        throw std::invalid_argument(
+            "mctimeseries ctor: expected 1-D array");
+    std::size_t n = static_cast<std::size_t>(view.shape(0));
+    std::vector<ValueType> out(n);
+    double const * data = view.data();
+    for (std::size_t i = 0; i < n; ++i)
+        out[i] = static_cast<ValueType>(data[i]);
+    return out;
+}
+// __repr__ helper for any ALPS type with an ostream operator.
+template <typename T>
+std::string stream_repr(T const & x) {
+    std::ostringstream ss;
+    ss << x;
+    return ss.str();
+}
+template <typename T>
+std::string value_with_error_repr(alps::alea::value_with_error<T> const & v) {
+    std::ostringstream ss;
+    ss << v.mean() << " +/- " << v.error();
+    return ss.str();
+}
+// Numpy-returning wrappers for vector-valued alps::alea free
+// functions (mean, variance, uncorrelated_error, binning_error) —
+// the scalar overloads are bound directly and nanobind casts
+// their `double` return to a Python float automatically.
+template <typename T>
+nb::object mean_vector(T const & x) {
+    return seq_to_numpy(alps::alea::mean(x));
+}
+template <typename T>
+nb::object variance_vector(T const & x) {
+    return seq_to_numpy(alps::alea::variance(x));
+}
+// mctimeseries.timeseries() returns std::vector<ValueType>; hand
+// back to Python as numpy. For scalar ValueType we pack 1-D; for
+// vector<double> ValueType we pack 2-D. mctimeseries_view has the
+// same surface.
+template <typename ValueType, typename TS>
+nb::object ts_to_numpy_scalar(TS const & ts) {
+    auto const & v = ts.timeseries();
+    return seq_to_numpy(v);
+}
+template <typename TS>
+nb::object ts_to_numpy_vector_rows(TS const & ts) {
+    auto const & rows = ts.timeseries();
+    if (rows.empty())
+        return alps::python::make_numpy_array<double>(
+            static_cast<double const*>(nullptr), {std::size_t{0}, std::size_t{0}});
+    std::size_t nrows = rows.size();
+    std::size_t ncols = rows.front().size();
+    for (auto const & row : rows)
+        if (row.size() != ncols)
+            throw std::runtime_error("mctimeseries has ragged rows; cannot shape as numpy 2-D");
+    std::vector<double> flat(nrows * ncols);
+    double * dst = flat.data();
+    for (auto const & row : rows) {
+        for (std::size_t j = 0; j < ncols; ++j)
+            *dst++ = static_cast<double>(row[j]);
+    }
+    return alps::python::make_numpy_array<double>(flat.data(), {nrows, ncols});
+}
+} // namespace
+NB_MODULE(pyalea_c, m) {
+    m.doc() = "ALPS alea bindings (nanobind)";
+    // ─── scalar-valarray observables ─────────────────────────────────
+    using RealVecObs = alps::alea::WrappedValarrayObservable<alps::RealVectorObservable>;
+    using RealVecTsObs = alps::alea::WrappedValarrayObservable<alps::RealVectorTimeSeriesObservable>;
+    #define ALPS_PY_EXPORT_VECTOROBSERVABLE(Wrapper, PyName)                              \
+        nb::class_<Wrapper>(m, PyName)                                                    \
+            .def("__init__",                                                              \
+                 [](Wrapper * self, std::string name, int bins) {                         \
+                     new (self) Wrapper(name, bins);                                      \
+                 },                                                                       \
+                 nb::arg("name"), nb::arg("bins") = 0)                                    \
+            .def("__repr__",     &Wrapper::representation)                                \
+            .def("__deepcopy__",                                                          \
+                 [](Wrapper const & self, nb::handle /*memo*/) {                          \
+                     return Wrapper(self);                                                \
+                 })                                                                       \
+            .def("__lshift__",   &Wrapper::push, nb::arg("array"))                        \
+            .def("save",         &Wrapper::save, nb::arg("filename"))                     \
+            .def_prop_ro("mean",     &Wrapper::mean)                                      \
+            .def_prop_ro("error",    &Wrapper::error)                                     \
+            .def_prop_ro("tau",      &Wrapper::tau)                                       \
+            .def_prop_ro("variance", &Wrapper::variance)                                  \
+            .def_prop_ro("count",    &Wrapper::count)                                     \
+            .def_prop_ro("converged_errors", &Wrapper::converged_errors)
+    ALPS_PY_EXPORT_VECTOROBSERVABLE(RealVecObs,   "RealVectorObservable");
+    ALPS_PY_EXPORT_VECTOROBSERVABLE(RealVecTsObs, "RealVectorTimeSeriesObservable");
+    #undef ALPS_PY_EXPORT_VECTOROBSERVABLE
+    // ─── scalar simple observables ───────────────────────────────────
+    #define ALPS_PY_EXPORT_SIMPLEOBSERVABLE(AlpsClass, PyName)                            \
+        nb::class_<alps::AlpsClass>(m, PyName)                                            \
+            .def("__init__",                                                              \
+                 [](alps::AlpsClass * self, std::string name, int bins) {                 \
+                     new (self) alps::AlpsClass(name, bins);                              \
+                 },                                                                       \
+                 nb::arg("name"), nb::arg("bins") = 0)                                    \
+            .def("__deepcopy__",                                                          \
+                 [](alps::AlpsClass const & self, nb::handle /*memo*/) {                  \
+                     return alps::AlpsClass(self);                                        \
+                 })                                                                       \
+            .def("__repr__",   &alps::AlpsClass::representation)                          \
+            .def("__lshift__", &alps::AlpsClass::operator<<)                              \
+            .def("save",       &alps::python::save_observable_to_hdf5<alps::AlpsClass>,   \
+                 nb::arg("filename"))                                                     \
+            .def_prop_ro("mean", &alps::AlpsClass::mean)                                  \
+            .def_prop_ro("error",                                                         \
+                static_cast<alps::AlpsClass::result_type (alps::AlpsClass::*)() const>(   \
+                    &alps::AlpsClass::error))                                             \
+            .def_prop_ro("tau",      &alps::AlpsClass::tau)                               \
+            .def_prop_ro("variance", &alps::AlpsClass::variance)                          \
+            .def_prop_ro("count",    &alps::AlpsClass::count)                             \
+            .def_prop_ro("converged_errors", &alps::AlpsClass::converged_errors)
+    ALPS_PY_EXPORT_SIMPLEOBSERVABLE(RealObservable,           "RealObservable");
+    ALPS_PY_EXPORT_SIMPLEOBSERVABLE(RealTimeSeriesObservable, "RealTimeSeriesObservable");
+    #undef ALPS_PY_EXPORT_SIMPLEOBSERVABLE
+    // ─── value_with_error ────────────────────────────────────────────
+    nb::class_<alps::alea::value_with_error<double>>(m, "ValueWithError")
+        .def(nb::init<double, double>(),
+             nb::arg("mean") = 0.0, nb::arg("error") = 0.0)
+        .def_prop_ro("mean",  &alps::alea::value_with_error<double>::mean)
+        .def_prop_ro("error", &alps::alea::value_with_error<double>::error)
+        .def("__repr__", &value_with_error_repr<double>);
+    // ─── StdPairDouble ─────────────────────────────────────────────
+    // nanobind's STL caster already registered std::pair<double, double> as a
+    // Python tuple converter via the stl.h header; nanobind takes the
+    // same path via <nanobind/stl/pair.h>. Binding it again as a
+    // class_ would fight that, so use a thin attribute-access
+    // wrapper instead, and give integrated_autocorrelation_time a
+    // Python-side signature that accepts either StdPairDouble or a
+    // plain (float, float) tuple.
+    struct StdPairDouble {
+        double first{0.0};
+        double second{0.0};
+        StdPairDouble() = default;
+        StdPairDouble(double f, double s) : first(f), second(s) {}
+        operator std::pair<double, double>() const { return {first, second}; }
+    };
+    nb::class_<StdPairDouble>(m, "StdPairDouble",
+        "Pair of (fit slope, fit intercept) returned by the autocorrelation fit helpers.")
+        .def(nb::init<>())
+        .def(nb::init<double, double>(), nb::arg("first"), nb::arg("second"))
+        .def_rw("first",  &StdPairDouble::first)
+        .def_rw("second", &StdPairDouble::second)
+        .def("__repr__", [](StdPairDouble const & p) {
+            std::ostringstream ss;
+            ss << "StdPairDouble(" << p.first << ", " << p.second << ")";
+            return ss.str();
+        });
+    // ─── mctimeseries<T> / mctimeseries_view<T> bindings ────────────
+    //
+    // Numpy-ctor: take any handle and copy through numpy_to_vector.
+    // The timeseries() method returns numpy arrays directly; the
+    // 2-D overload for vector<double> goes through the
+    // ts_to_numpy_vector_rows helper.
+    #define ALPS_PY_EXPORT_MCTIMESERIES_SCALAR(Value, PyName)                             \
+        nb::class_<alps::alea::mctimeseries<Value>>(m, PyName)                            \
+            .def(nb::init<>())                                                            \
+            .def("__init__",                                                              \
+                 [](alps::alea::mctimeseries<Value> * self, nb::handle a) {               \
+                     new (self) alps::alea::mctimeseries<Value>(                          \
+                         numpy_to_vector<Value>(a));                                      \
+                 })                                                                       \
+            .def(nb::init<alps::alea::mcdata<Value>>())                                   \
+            .def("timeseries", [](alps::alea::mctimeseries<Value> const & self) {         \
+                    return ts_to_numpy_scalar<Value>(self);                               \
+                })                                                                        \
+            .def_prop_ro("size", &alps::alea::mctimeseries<Value>::size)                  \
+            .def("__repr__", &stream_repr<alps::alea::mctimeseries<Value>>);              \
+        nb::class_<alps::alea::mctimeseries_view<Value>>(m, PyName "View")                \
+            .def(nb::init<alps::alea::mctimeseries<Value>>())                             \
+            .def(nb::init<alps::alea::mctimeseries_view<Value>>())                        \
+            .def("timeseries", [](alps::alea::mctimeseries_view<Value> const & self) {    \
+                    return ts_to_numpy_scalar<Value>(self);                               \
+                })                                                                        \
+            .def_prop_ro("size", &alps::alea::mctimeseries_view<Value>::size)             \
+            .def("__repr__", &stream_repr<alps::alea::mctimeseries_view<Value>>)
+    ALPS_PY_EXPORT_MCTIMESERIES_SCALAR(double, "MCScalarTimeseries");
+    #undef ALPS_PY_EXPORT_MCTIMESERIES_SCALAR
+    // Vector-valued mctimeseries: ctor from a 2-D numpy array, rows
+    // are time samples. timeseries() returns 2-D.
+    using VecTs    = alps::alea::mctimeseries<std::vector<double>>;
+    using VecTsV   = alps::alea::mctimeseries_view<std::vector<double>>;
+    using VecMcD   = alps::alea::mcdata<std::vector<double>>;
+    nb::class_<VecTs>(m, "MCVectorTimeseries")
+        .def(nb::init<>())
+        .def("__init__", [](VecTs * self, nb::handle a) {
+                auto view = alps::python::as_contiguous<double>(a);
+                if (view.ndim() != 2)
+                    throw std::invalid_argument(
+                        "MCVectorTimeseries ctor: expected 2-D array");
+                std::size_t nrows = static_cast<std::size_t>(view.shape(0));
+                std::size_t ncols = static_cast<std::size_t>(view.shape(1));
+                std::vector<std::vector<double>> rows(nrows);
+                double const * data = view.data();
+                for (std::size_t i = 0; i < nrows; ++i) {
+                    rows[i].resize(ncols);
+                    for (std::size_t j = 0; j < ncols; ++j)
+                        rows[i][j] = data[i * ncols + j];
+                }
+                new (self) VecTs(rows);
+            })
+        .def(nb::init<VecMcD>())
+        .def("timeseries", [](VecTs const & self) { return ts_to_numpy_vector_rows(self); })
+        .def_prop_ro("size", &VecTs::size)
+        .def("__repr__", &stream_repr<VecTs>);
+    nb::class_<VecTsV>(m, "MCVectorTimeseriesView")
+        .def(nb::init<VecTs>())
+        .def(nb::init<VecTsV>())
+        .def("timeseries", [](VecTsV const & self) { return ts_to_numpy_vector_rows(self); })
+        .def_prop_ro("size", &VecTsV::size)
+        .def("__repr__", &stream_repr<VecTsV>);
+    // ─── alps::alea free functions over mcdata / mctimeseries ────────
+    #define DEF_ALL(name, fn)                                                             \
+        /* scalar-valued */                                                               \
+        m.def(name, static_cast<double (*)(alps::alea::mcdata<double> const &)>(&fn));    \
+        m.def(name, static_cast<double (*)(alps::alea::mctimeseries<double> const &)>(&fn)); \
+        m.def(name, static_cast<double (*)(alps::alea::mctimeseries_view<double> const &)>(&fn));
+    // size — works for both scalar and vector value types.
+    m.def("size", static_cast<std::size_t (*)(alps::alea::mcdata<double> const &)>(&alps::size));
+    m.def("size", static_cast<std::size_t (*)(alps::alea::mctimeseries<double> const &)>(&alps::size));
+    m.def("size", static_cast<std::size_t (*)(alps::alea::mctimeseries_view<double> const &)>(&alps::size));
+    m.def("size", static_cast<std::size_t (*)(alps::alea::mcdata<std::vector<double>> const &)>(&alps::size));
+    m.def("size", static_cast<std::size_t (*)(alps::alea::mctimeseries<std::vector<double>> const &)>(&alps::size));
+    m.def("size", static_cast<std::size_t (*)(alps::alea::mctimeseries_view<std::vector<double>> const &)>(&alps::size));
+    // mean — scalar overloads return double, vector overloads return numpy.
+    DEF_ALL("mean", alps::alea::mean)
+    m.def("mean", &mean_vector<alps::alea::mcdata<std::vector<double>>>);
+    m.def("mean", &mean_vector<alps::alea::mctimeseries<std::vector<double>>>);
+    m.def("mean", &mean_vector<alps::alea::mctimeseries_view<std::vector<double>>>);
+    // variance — same pattern.
+    DEF_ALL("variance", alps::alea::variance)
+    m.def("variance", &variance_vector<alps::alea::mcdata<std::vector<double>>>);
+    m.def("variance", &variance_vector<alps::alea::mctimeseries<std::vector<double>>>);
+    m.def("variance", &variance_vector<alps::alea::mctimeseries_view<std::vector<double>>>);
+    // integrated_autocorrelation_time — scalar only. The C++ signature
+    // takes the (slope, intercept) pair by const-ref.
+    m.def("integrated_autocorrelation_time",
+          static_cast<double (*)(alps::alea::mctimeseries<double> const &,
+                                 std::pair<double, double> const &)>(
+              &alps::alea::integrated_autocorrelation_time));
+    m.def("integrated_autocorrelation_time",
+          static_cast<double (*)(alps::alea::mctimeseries_view<double> const &,
+                                 std::pair<double, double> const &)>(
+              &alps::alea::integrated_autocorrelation_time));
+    // running_mean / reverse_running_mean — scalar only, mctimeseries-valued.
+    m.def("running_mean",
+          static_cast<alps::alea::mctimeseries<double> (*)(alps::alea::mcdata<double> const &)>(
+              &alps::alea::running_mean));
+    m.def("running_mean",
+          static_cast<alps::alea::mctimeseries<double> (*)(alps::alea::mctimeseries<double> const &)>(
+              &alps::alea::running_mean));
+    m.def("running_mean",
+          static_cast<alps::alea::mctimeseries<double> (*)(alps::alea::mctimeseries_view<double> const &)>(
+              &alps::alea::running_mean));
+    m.def("reverse_running_mean",
+          static_cast<alps::alea::mctimeseries<double> (*)(alps::alea::mcdata<double> const &)>(
+              &alps::alea::reverse_running_mean));
+    m.def("reverse_running_mean",
+          static_cast<alps::alea::mctimeseries<double> (*)(alps::alea::mctimeseries<double> const &)>(
+              &alps::alea::reverse_running_mean));
+    m.def("reverse_running_mean",
+          static_cast<alps::alea::mctimeseries<double> (*)(alps::alea::mctimeseries_view<double> const &)>(
+              &alps::alea::reverse_running_mean));
+    #undef DEF_ALL
+}
diff --git a/bindings/python/pyalps/cpp/pymcdata.cpp b/bindings/python/pyalps/cpp/pymcdata.cpp
new file mode 100644
index 000000000..017299470
--- /dev/null
+++ b/bindings/python/pyalps/cpp/pymcdata.cpp
@@ -0,0 +1,329 @@
+// Copyright (C) 1994-2010 by Ping Nang Ma <pingnang@itp.phys.ethz.ch>,
+//               Lukas Gamper <gamperl@gmail.com>,
+//               Matthias Troyer <troyer@itp.phys.ethz.ch>
+//               2026       by the ALPS collaboration
+// Part of the ALPS Project — see LICENSE.txt for full license text.
+// SPDX-License-Identifier: MIT
+#include <nanobind/nanobind.h>
+#include <nanobind/stl/string.h>
+#include <nanobind/stl/vector.h>
+#include <alps/alea/mcdata.hpp>
+#include "numpy_compat.hpp"
+#include <algorithm>
+#include <cmath>
+#include <cstddef>
+#include <sstream>
+#include <string>
+#include <utility>
+#include <vector>
+namespace nb = nanobind;
+namespace alps {
+    namespace python {
+        // Build a 1-D numpy.ndarray (dtype=float64) from a
+        // std::vector<double>. Allocates a fresh array via
+        // numpy.empty + memcpy through the buffer protocol — no
+        // numpy headers.
+        inline nb::object vec_to_numpy(std::vector<double> const & v) {
+            return make_numpy_array<double>(v.data(), {v.size()});
+        }
+        // Build a 2-D numpy.ndarray from a vector<vector<double>>.
+        // Rows must be equal-length; on mismatch throw a value error
+        // (mcdata doesn't produce ragged bins/jackknife tables).
+        inline nb::object matrix_to_numpy(std::vector<std::vector<double>> const & m) {
+            std::size_t rows = m.size();
+            std::size_t cols = rows > 0 ? m.front().size() : 0;
+            for (auto const & row : m) {
+                if (row.size() != cols)
+                    throw std::runtime_error("mcdata returned ragged 2-D table; refusing to convert to numpy");
+            }
+            std::vector<double> flat(rows * cols);
+            double * dst = flat.data();
+            for (auto const & row : m) {
+                std::copy(row.begin(), row.end(), dst);
+                dst += cols;
+            }
+            return make_numpy_array<double>(flat.data(), {rows, cols});
+        }
+        // __repr__ for mcdata<double> — "<mean> +/- <error>".
+        template <typename T>
+        std::string print_mcdata_scalar(alps::alea::mcdata<T> const & self) {
+            std::ostringstream ss;
+            ss << self.mean() << " +/- " << self.error();
+            return ss.str();
+        }
+        // __repr__ for mcdata<vector<T>> — newline-joined scalar reprs.
+        // The const_iterator's operator-> returns through boost::addressof
+        // on a rvalue scalar mcdata<T>, which triggers the deleted
+        // overload; take a local copy instead.
+        template <typename T>
+        std::string print_mcdata_vector(alps::alea::mcdata<std::vector<T>> const & self) {
+            std::ostringstream ss;
+            bool first = true;
+            for (auto it = self.begin(); it != self.end(); ++it) {
+                if (!first) ss << "\n";
+                first = false;
+                alps::alea::mcdata<T> entry = *it;
+                ss << entry.mean() << " +/- " << entry.error();
+            }
+            return ss.str();
+        }
+        // __format__ for mcdata<T> — defers to builtins.format on
+        // mean and error separately, then joins with " +/- ".
+        template <typename T>
+        std::string format_mcdata_scalar(alps::alea::mcdata<T> const & self,
+                                         std::string const & spec) {
+            nb::object fmt = nb::module_::import_("builtins").attr("format");
+            std::string m = nb::cast<std::string>(fmt(self.mean(), spec));
+            std::string e = nb::cast<std::string>(fmt(self.error(), spec));
+            return m + " +/- " + e;
+        }
+        template <typename T>
+        std::string format_mcdata_vector(alps::alea::mcdata<std::vector<T>> const & self,
+                                         std::string const & spec) {
+            std::ostringstream ss;
+            bool first = true;
+            for (auto it = self.begin(); it != self.end(); ++it) {
+                if (!first) ss << "\n";
+                first = false;
+                alps::alea::mcdata<T> entry = *it;
+                ss << format_mcdata_scalar(entry, spec);
+            }
+            return ss.str();
+        }
+        // Indexing: mcdata<vector<T>>[i] returns a scalar mcdata<T>
+        // view; mcdata<vector<T>>[slice] returns a sliced
+        // mcdata<vector<T>>.
+        template <typename T>
+        nb::object mcdata_vector_getitem(alps::alea::mcdata<std::vector<T>> & data,
+                                         nb::object const & key) {
+            std::size_t n = data.mean().size();
+            if (nb::isinstance<nb::slice>(key)) {
+                nb::slice s = nb::borrow<nb::slice>(key);
+                auto [start, stop, step, slicelength] = s.compute(n);
+                if (step != 1)
+                    throw nb::index_error("slice step size not supported.");
+                if (start > stop)
+                    return nb::cast(alps::alea::mcdata<std::vector<T>>());
+                return nb::cast(alps::alea::mcdata<std::vector<T>>(
+                    data, start, stop));
+            }
+            long index = nb::cast<long>(key);
+            if (index < 0)
+                index += static_cast<long>(n);
+            if (index < 0 || static_cast<std::size_t>(index) >= n)
+                throw nb::index_error("Index out of range");
+            return nb::cast(alps::alea::mcdata<T>(
+                data, static_cast<std::size_t>(index)));
+        }
+        template <typename T>
+        bool mcdata_vector_contains(alps::alea::mcdata<std::vector<T>> & data,
+                                    nb::object const & key) {
+            // Best-effort: accept either a scalar mcdata<T> or a value
+            // that casts cleanly; anything else is "not in".
+            // mcdata<std::vector<T>>::const_iterator does not fully
+            // model std::ranges::input_range (the iterator traits
+            // needed for borrowed_iterator_t aren't wired up), so the
+            // ranges-form of std::find is ill-formed here. Keep the
+            // classical iterator pair.
+            try {
+                auto probe = nb::cast<alps::alea::mcdata<T>>(key);
+                return std::find(data.begin(), data.end(), probe) != data.end();
+            } catch (nb::cast_error const &) {
+            }
+            return false;
+        }
+    }
+}
+NB_MODULE(pymcdata_c, m) {
+    using alps::alea::mcdata;
+    namespace pymod = alps::python;
+    using Scalar = mcdata<double>;
+    using Vector = mcdata<std::vector<double>>;
+    // mcdata's transcendentals (sq, cb, sqrt, cbrt, exp, log, sin, …)
+    // live in namespace alps::alea and are found via ADL when we pass
+    // an mcdata<T> argument. Capture each as a lambda so the binding
+    // doesn't have to cast through overloaded name lookup, which
+    // fights <cmath>'s own abs/pow/sqrt etc. sitting in global scope.
+    nb::class_<Scalar>(m, "MCScalarData",
+        "Scalar Monte Carlo data. Supports +, -, *, /, +=, -=, *=, /=, "
+        "abs, pow, sq, cb, sqrt, cbrt, exp, log, sin, cos, tan, sinh, cosh, tanh.")
+        .def(nb::init<>())
+        .def(nb::init<double>(), nb::arg("mean"))
+        .def(nb::init<double, double>(), nb::arg("mean"), nb::arg("error"))
+        .def_prop_ro("mean",     [](Scalar const & v) { return v.mean(); })
+        .def_prop_ro("error",    [](Scalar const & v) { return v.error(); })
+        .def_prop_ro("tau",      [](Scalar const & v) { return v.tau(); })
+        .def_prop_ro("variance", [](Scalar const & v) { return v.variance(); })
+        .def_prop_ro("bins",     [](Scalar const & v) {
+                                      return pymod::vec_to_numpy(v.bins());
+                                  })
+        .def_prop_ro("jackknife",[](Scalar const & v) {
+                                      return pymod::vec_to_numpy(v.jackknife());
+                                  })
+        .def_prop_ro("count",    &Scalar::count)
+        .def("__repr__",   &pymod::print_mcdata_scalar<double>)
+        .def("__format__", &pymod::format_mcdata_scalar<double>,
+             nb::arg("format_spec"))
+        .def("__deepcopy__",
+             [](Scalar const & self, nb::handle /*memo*/) {
+                 return Scalar(self);
+             })
+        .def("__abs__", [](Scalar x) { using alps::alea::abs; return abs(std::move(x)); })
+        .def("__pow__", [](Scalar x, double e) { using alps::alea::pow; return pow(std::move(x), e); })
+        // Unary - / + on mcdata produce new values; wrap manually
+        // because the library's operator+()/-() signatures aren't
+        // const-returning (which is what nb::self expects).
+        .def("__pos__", [](Scalar self) { return +self; })
+        .def("__neg__", [](Scalar self) { return -self; })
+        // In-place operators — modify self in place, return reference.
+        .def("__iadd__", [](Scalar & s, Scalar const & o) -> Scalar & { return s += o; }, nb::is_operator())
+        .def("__iadd__", [](Scalar & s, double o)         -> Scalar & { return s += o; }, nb::is_operator())
+        .def("__isub__", [](Scalar & s, Scalar const & o) -> Scalar & { return s -= o; }, nb::is_operator())
+        .def("__isub__", [](Scalar & s, double o)         -> Scalar & { return s -= o; }, nb::is_operator())
+        .def("__imul__", [](Scalar & s, Scalar const & o) -> Scalar & { return s *= o; }, nb::is_operator())
+        .def("__imul__", [](Scalar & s, double o)         -> Scalar & { return s *= o; }, nb::is_operator())
+        .def("__itruediv__", [](Scalar & s, Scalar const & o) -> Scalar & { return s /= o; }, nb::is_operator())
+        .def("__itruediv__", [](Scalar & s, double o)         -> Scalar & { return s /= o; }, nb::is_operator())
+        // Binary +/-/*//: forward and reflected forms. nb::is_operator()
+        // marks them so mixed-type failures return NotImplemented and
+        // Python's reflected operator machinery gets a turn.
+        .def("__add__",  [](Scalar const & a, Scalar const & b) { return a + b; }, nb::is_operator())
+        .def("__radd__", [](Scalar const & a, Scalar const & b) { return b + a; }, nb::is_operator())
+        .def("__add__",  [](Scalar const & a, double b)         { return a + b; }, nb::is_operator())
+        .def("__radd__", [](Scalar const & a, double b)         { return b + a; }, nb::is_operator())
+        .def("__sub__",  [](Scalar const & a, Scalar const & b) { return a - b; }, nb::is_operator())
+        .def("__rsub__", [](Scalar const & a, Scalar const & b) { return b - a; }, nb::is_operator())
+        .def("__sub__",  [](Scalar const & a, double b)         { return a - b; }, nb::is_operator())
+        .def("__rsub__", [](Scalar const & a, double b)         { return b - a; }, nb::is_operator())
+        .def("__mul__",  [](Scalar const & a, Scalar const & b) { return a * b; }, nb::is_operator())
+        .def("__rmul__", [](Scalar const & a, Scalar const & b) { return b * a; }, nb::is_operator())
+        .def("__mul__",  [](Scalar const & a, Vector const & b) { return a * b; }, nb::is_operator())
+        .def("__rmul__", [](Scalar const & a, Vector const & b) { return b * a; }, nb::is_operator())
+        .def("__mul__",  [](Scalar const & a, double b)         { return a * b; }, nb::is_operator())
+        .def("__rmul__", [](Scalar const & a, double b)         { return b * a; }, nb::is_operator())
+        .def("__truediv__",  [](Scalar const & a, Scalar const & b) { return a / b; }, nb::is_operator())
+        .def("__rtruediv__", [](Scalar const & a, Scalar const & b) { return b / a; }, nb::is_operator())
+        .def("__truediv__",  [](Scalar const & a, double b)         { return a / b; }, nb::is_operator())
+        .def("__rtruediv__", [](Scalar const & a, double b)         { return b / a; }, nb::is_operator())
+        .def("sq",   [](Scalar x) { using alps::alea::sq;   return sq(std::move(x));   })
+        .def("cb",   [](Scalar x) { using alps::alea::cb;   return cb(std::move(x));   })
+        .def("sqrt", [](Scalar x) { using alps::alea::sqrt; return sqrt(std::move(x)); })
+        .def("cbrt", [](Scalar x) { using alps::alea::cbrt; return cbrt(std::move(x)); })
+        .def("exp",  [](Scalar x) { using alps::alea::exp;  return exp(std::move(x));  })
+        .def("log",  [](Scalar x) { using alps::alea::log;  return log(std::move(x));  })
+        .def("sin",  [](Scalar x) { using alps::alea::sin;  return sin(std::move(x));  })
+        .def("cos",  [](Scalar x) { using alps::alea::cos;  return cos(std::move(x));  })
+        .def("tan",  [](Scalar x) { using alps::alea::tan;  return tan(std::move(x));  })
+        .def("sinh", [](Scalar x) { using alps::alea::sinh; return sinh(std::move(x)); })
+        .def("cosh", [](Scalar x) { using alps::alea::cosh; return cosh(std::move(x)); })
+        .def("tanh", [](Scalar x) { using alps::alea::tanh; return tanh(std::move(x)); })
+        .def("set_bin_size",   &Scalar::set_bin_size)
+        .def("set_bin_number", &Scalar::set_bin_number)
+        .def("discard_bins",   &Scalar::discard_bins)
+        .def("merge", static_cast<void (Scalar::*)(Scalar const &)>(&Scalar::merge))
+        .def("save", static_cast<void (Scalar::*)(std::string const &, std::string const &) const>(&Scalar::save),
+             nb::arg("filename"), nb::arg("observable_name"))
+        .def("load", static_cast<void (Scalar::*)(std::string const &, std::string const &)>(&Scalar::load),
+             nb::arg("filename"), nb::arg("observable_name"));
+    nb::class_<Vector>(m, "MCVectorData",
+        "Vector-valued Monte Carlo data.")
+        .def(nb::init<>())
+        .def(nb::init<std::vector<double>>(), nb::arg("mean"))
+        .def(nb::init<std::vector<double>, std::vector<double>>(),
+             nb::arg("mean"), nb::arg("error"))
+        .def("__len__",
+             [](Vector & v) {
+                 return v.mean().size();
+             })
+        .def("__getitem__", &pymod::mcdata_vector_getitem<double>)
+        .def("__contains__", &pymod::mcdata_vector_contains<double>)
+        .def_prop_ro("mean",     [](Vector const & v) {
+                                      return pymod::vec_to_numpy(v.mean());
+                                  })
+        .def_prop_ro("error",    [](Vector const & v) {
+                                      return pymod::vec_to_numpy(v.error());
+                                  })
+        .def_prop_ro("tau",      [](Vector const & v) {
+                                      return pymod::vec_to_numpy(v.tau());
+                                  })
+        .def_prop_ro("variance", [](Vector const & v) {
+                                      return pymod::vec_to_numpy(v.variance());
+                                  })
+        .def_prop_ro("bins",     [](Vector const & v) {
+                                      return pymod::matrix_to_numpy(v.bins());
+                                  })
+        .def_prop_ro("jackknife",[](Vector const & v) {
+                                      return pymod::matrix_to_numpy(v.jackknife());
+                                  })
+        .def_prop_ro("count", &Vector::count)
+        .def("__repr__",   &pymod::print_mcdata_vector<double>)
+        .def("__format__", &pymod::format_mcdata_vector<double>,
+             nb::arg("format_spec"))
+        .def("__deepcopy__",
+             [](Vector const & self, nb::handle /*memo*/) {
+                 return Vector(self);
+             })
+        .def("__abs__", [](Vector x) { using alps::alea::abs; return abs(std::move(x)); })
+        .def("__pow__", [](Vector x, double e) { using alps::alea::pow; return pow(std::move(x), e); })
+        .def("__pos__", [](Vector self) { return +self; })
+        .def("__neg__", [](Vector self) { return -self; })
+        .def("__eq__",  [](Vector const & a, Vector const & b) { return a == b; }, nb::is_operator())
+        // In-place operators.
+        .def("__iadd__", [](Vector & s, Vector const & o)              -> Vector & { return s += o; }, nb::is_operator())
+        .def("__iadd__", [](Vector & s, std::vector<double> const & o) -> Vector & { return s += o; }, nb::is_operator())
+        .def("__isub__", [](Vector & s, Vector const & o)              -> Vector & { return s -= o; }, nb::is_operator())
+        .def("__isub__", [](Vector & s, std::vector<double> const & o) -> Vector & { return s -= o; }, nb::is_operator())
+        .def("__imul__", [](Vector & s, Vector const & o)              -> Vector & { return s *= o; }, nb::is_operator())
+        .def("__imul__", [](Vector & s, std::vector<double> const & o) -> Vector & { return s *= o; }, nb::is_operator())
+        .def("__itruediv__", [](Vector & s, Vector const & o)              -> Vector & { return s /= o; }, nb::is_operator())
+        .def("__itruediv__", [](Vector & s, std::vector<double> const & o) -> Vector & { return s /= o; }, nb::is_operator())
+        // Binary operators — Vector ↔ Vector / Scalar / vector<double> / double.
+        .def("__add__",  [](Vector const & a, Vector const & b)              { return a + b; }, nb::is_operator())
+        .def("__radd__", [](Vector const & a, Vector const & b)              { return b + a; }, nb::is_operator())
+        .def("__add__",  [](Vector const & a, std::vector<double> const & b) { return a + b; }, nb::is_operator())
+        .def("__radd__", [](Vector const & a, std::vector<double> const & b) { return b + a; }, nb::is_operator())
+        .def("__sub__",  [](Vector const & a, Vector const & b)              { return a - b; }, nb::is_operator())
+        .def("__rsub__", [](Vector const & a, Vector const & b)              { return b - a; }, nb::is_operator())
+        .def("__sub__",  [](Vector const & a, std::vector<double> const & b) { return a - b; }, nb::is_operator())
+        .def("__rsub__", [](Vector const & a, std::vector<double> const & b) { return b - a; }, nb::is_operator())
+        .def("__mul__",  [](Vector const & a, Vector const & b)              { return a * b; }, nb::is_operator())
+        .def("__mul__",  [](Vector const & a, Scalar const & b)              { return a * b; }, nb::is_operator())
+        .def("__rmul__", [](Vector const & a, Scalar const & b)              { return b * a; }, nb::is_operator())
+        .def("__rmul__", [](Vector const & a, Vector const & b)              { return b * a; }, nb::is_operator())
+        .def("__mul__",  [](Vector const & a, std::vector<double> const & b) { return a * b; }, nb::is_operator())
+        .def("__rmul__", [](Vector const & a, std::vector<double> const & b) { return b * a; }, nb::is_operator())
+        .def("__truediv__",  [](Vector const & a, Vector const & b)              { return a / b; }, nb::is_operator())
+        .def("__truediv__",  [](Vector const & a, Scalar const & b)              { return a / b; }, nb::is_operator())
+        .def("__rtruediv__", [](Vector const & a, Vector const & b)              { return b / a; }, nb::is_operator())
+        .def("__truediv__",  [](Vector const & a, std::vector<double> const & b) { return a / b; }, nb::is_operator())
+        .def("__rtruediv__", [](Vector const & a, std::vector<double> const & b) { return b / a; }, nb::is_operator())
+        // Vector ↔ double.
+        .def("__add__",  [](Vector const & a, double b) { return a + b; }, nb::is_operator())
+        .def("__radd__", [](Vector const & a, double b) { return b + a; }, nb::is_operator())
+        .def("__sub__",  [](Vector const & a, double b) { return a - b; }, nb::is_operator())
+        .def("__rsub__", [](Vector const & a, double b) { return b - a; }, nb::is_operator())
+        .def("__mul__",  [](Vector const & a, double b) { return a * b; }, nb::is_operator())
+        .def("__rmul__", [](Vector const & a, double b) { return b * a; }, nb::is_operator())
+        .def("__truediv__",  [](Vector const & a, double b) { return a / b; }, nb::is_operator())
+        .def("__rtruediv__", [](Vector const & a, double b) { return b / a; }, nb::is_operator())
+        .def("sq",   [](Vector x) { using alps::alea::sq;   return sq(std::move(x));   })
+        .def("cb",   [](Vector x) { using alps::alea::cb;   return cb(std::move(x));   })
+        .def("sqrt", [](Vector x) { using alps::alea::sqrt; return sqrt(std::move(x)); })
+        .def("cbrt", [](Vector x) { using alps::alea::cbrt; return cbrt(std::move(x)); })
+        .def("exp",  [](Vector x) { using alps::alea::exp;  return exp(std::move(x));  })
+        .def("log",  [](Vector x) { using alps::alea::log;  return log(std::move(x));  })
+        .def("sin",  [](Vector x) { using alps::alea::sin;  return sin(std::move(x));  })
+        .def("cos",  [](Vector x) { using alps::alea::cos;  return cos(std::move(x));  })
+        .def("tan",  [](Vector x) { using alps::alea::tan;  return tan(std::move(x));  })
+        .def("sinh", [](Vector x) { using alps::alea::sinh; return sinh(std::move(x)); })
+        .def("cosh", [](Vector x) { using alps::alea::cosh; return cosh(std::move(x)); })
+        .def("tanh", [](Vector x) { using alps::alea::tanh; return tanh(std::move(x)); })
+        .def("set_bin_size",   &Vector::set_bin_size)
+        .def("set_bin_number", &Vector::set_bin_number)
+        .def("discard_bins",   &Vector::discard_bins)
+        .def("merge", static_cast<void (Vector::*)(Vector const &)>(&Vector::merge))
+        .def("save", static_cast<void (Vector::*)(std::string const &, std::string const &) const>(&Vector::save),
+             nb::arg("filename"), nb::arg("observable_name"))
+        .def("load", static_cast<void (Vector::*)(std::string const &, std::string const &)>(&Vector::load),
+             nb::arg("filename"), nb::arg("observable_name"));
+}
diff --git a/bindings/python/pyalps/cpp/pytools.cpp b/bindings/python/pyalps/cpp/pytools.cpp
new file mode 100644
index 000000000..3edc408ff
--- /dev/null
+++ b/bindings/python/pyalps/cpp/pytools.cpp
@@ -0,0 +1,52 @@
+// Copyright (C) 1994-2009 by Ping Nang Ma <pingnang@itp.phys.ethz.ch>,
+//               Matthias Troyer <troyer@itp.phys.ethz.ch>,
+//               Bela Bauer <bauerb@itp.phys.ethz.ch>
+//               2026       by the ALPS collaboration
+// Part of the ALPS Project — see LICENSE.txt for full license text.
+// SPDX-License-Identifier: MIT
+#include <nanobind/nanobind.h>
+#include <nanobind/stl/string.h>
+#include <alps/scheduler/convert.h>
+#include <alps/utility/encode.hpp>
+#include <alps/random.h>
+#include <alps/parser/xslt_path.h>
+#include <alps/random/mersenne_twister.hpp>
+#include <boost/random/uniform_01.hpp>
+#include <boost/random/variate_generator.hpp>
+namespace nb = nanobind;
+typedef boost::variate_generator<boost::mt19937, boost::uniform_01<double> > random_01;
+class WrappedRNG : public random_01
+{
+public:
+    WrappedRNG(int seed = 0)
+        : random_01(boost::mt19937(seed), boost::uniform_01<double>())
+    {
+    }
+};
+NB_MODULE(pytools_c, m) {
+    m.doc() = "ALPS tools bindings (nanobind)";
+    m.def("convert2xml",
+          &alps::convert2xml,
+          "Convert an ALPS file to XML. Returns the path to the XML file.");
+    m.def("hdf5_name_encode",
+          &alps::hdf5_name_encode,
+          "Escape a string for use inside an HDF5 path name.");
+    m.def("hdf5_name_decode",
+          &alps::hdf5_name_decode,
+          "Un-escape a string taken from an HDF5 path name.");
+    m.def("search_xml_library_path",
+          &alps::search_xml_library_path,
+          "Resolve an ALPS library XML / XSL file to its full path.");
+    nb::class_<WrappedRNG>(m, "rng",
+                           "Mersenne-Twister uniform random number generator in [0, 1).")
+        .def(nb::init<int>(), nb::arg("seed") = 0)
+        .def("__deepcopy__",
+             [](WrappedRNG const & self, nb::handle /*memo*/) {
+                 return WrappedRNG(self);
+             },
+             "Return a fresh copy of the RNG carrying the same state.")
+        .def("__call__",
+             static_cast<WrappedRNG::result_type (WrappedRNG::*)()>(
+                 &WrappedRNG::operator()),
+             "Return a uniform random number in [0, 1).");
+}
diff --git a/bindings/python/pyalps/cpp/save_observable_to_hdf5.hpp b/bindings/python/pyalps/cpp/save_observable_to_hdf5.hpp
new file mode 100644
index 000000000..556634686
--- /dev/null
+++ b/bindings/python/pyalps/cpp/save_observable_to_hdf5.hpp
@@ -0,0 +1,15 @@
+// Copyright (C) 2010 by Matthias Troyer <troyer@comp-phys.org>,
+// Part of the ALPS Project — see LICENSE.txt for full license text.
+// SPDX-License-Identifier: MIT
+#ifndef ALPS_PYTHON_VERY_LONG_FILENAME_FOR_SAVE_OBSERVABLE_TO_HDF5_HPP
+#define ALPS_PYTHON_VERY_LONG_FILENAME_FOR_SAVE_OBSERVABLE_TO_HDF5_HPP
+#include <alps/hdf5.hpp>
+namespace alps { namespace python {
+
+    template <typename Obs> void save_observable_to_hdf5(Obs const & obs, std::string const & filename) {
+        hdf5::archive ar(filename, "a");
+        ar["/simulation/results/"+obs.representation()] << obs;
+    }
+
+} } // end namespace alps::python
+#endif // ALPS_PYTHON_VERY_LONG_FILENAME_FOR_SAVE_OBSERVABLE_TO_HDF5_HPP
diff --git a/bindings/python/pyalps/pyproject.toml b/bindings/python/pyalps/pyproject.toml
new file mode 100644
index 000000000..1becb4ccb
--- /dev/null
+++ b/bindings/python/pyalps/pyproject.toml
@@ -0,0 +1,37 @@
+[build-system]
+requires = ["scikit-build-core>=1.0", "nanobind>=2.10"]
+build-backend = "scikit_build_core.build"
+
+[project]
+name = "pyalps"
+version = "2.3.4b1"
+description = "Python Applications and Libraries for Physics Simulations"
+readme = "README.md"
+requires-python = ">=3.10"
+license = "MIT"
+dependencies = ["numpy>=1.26", "scipy>=1.13"]
+
+[project.optional-dependencies]
+plot = ["matplotlib>=3.8"]
+test = ["pytest>=8"]
+
+[tool.scikit-build]
+cmake.source-dir = "."
+wheel.packages = ["src/pyalps"]
+wheel.force-include = { "LICENSE.txt" = "${SKBUILD_METADATA_DIR}/licenses/LICENSE.txt" }
+build.verbose = true
+
+[tool.scikit-build.cmake.define]
+ALPS_DIR = { env = "ALPS_DIR" }
+
+# Application bindings compile selected legacy application sources. Preserve
+# those sources when this subproject is distributed independently of the
+# repository checkout so wheels can also be rebuilt from the sdist.
+[tool.scikit-build.sdist.force-include]
+"../../../LICENSE.txt" = "LICENSE.txt"
+"../../../applications/dmft/qmc" = "_vendor/applications/dmft/qmc"
+"../../../applications/qmc/dwa" = "_vendor/applications/qmc/dwa"
+"../../../tool" = "_vendor/tool"
+
+[tool.cibuildwheel]
+manylinux-x86_64-image = "manylinux_2_28"
diff --git a/lib/pyalps/__init__.py b/bindings/python/pyalps/src/pyalps/__init__.py
similarity index 85%
rename from lib/pyalps/__init__.py
rename to bindings/python/pyalps/src/pyalps/__init__.py
index a88834ea8..c70dd433c 100644
--- a/lib/pyalps/__init__.py
+++ b/bindings/python/pyalps/src/pyalps/__init__.py
@@ -28,13 +28,6 @@
 # ****************************************************************************
 
 import sys
-import os.path
-if sys.platform == 'darwin' and not os.path.exists(os.path.expanduser('~/.matplotlib/matplotlibrc')):
-    try:
-        import matplotlib
-        matplotlib.use('macosx')
-    except ImportError:
-        pass
 
 from .dataset import *
 from .tools import *
@@ -42,5 +35,14 @@
 from .floatwitherror import FloatWithError
 from . import fit_wrapper
 
+# Optional solver modules are present when the wheel was built from an ALPS
+# checkout with application bindings enabled.
+try:
+    from ._ext import cthyb, ctint
+    sys.modules[__name__ + ".cthyb"] = cthyb
+    sys.modules[__name__ + ".ctint"] = ctint
+except ImportError:
+    pass
+
 # For ALPS DWA Application
 # from dwa import *
diff --git a/bindings/python/pyalps/src/pyalps/_ext/__init__.py b/bindings/python/pyalps/src/pyalps/_ext/__init__.py
new file mode 100644
index 000000000..f3bf869bc
--- /dev/null
+++ b/bindings/python/pyalps/src/pyalps/_ext/__init__.py
@@ -0,0 +1,4 @@
+# Copyright (C) 2026 by the ALPS collaboration
+# SPDX-License-Identifier: MIT
+
+"""Compiled nanobind extensions for pyalps."""
diff --git a/lib/pyalps/alea.py b/bindings/python/pyalps/src/pyalps/alea.py
similarity index 100%
rename from lib/pyalps/alea.py
rename to bindings/python/pyalps/src/pyalps/alea.py
diff --git a/lib/pyalps/alea_detail.py b/bindings/python/pyalps/src/pyalps/alea_detail.py
similarity index 100%
rename from lib/pyalps/alea_detail.py
rename to bindings/python/pyalps/src/pyalps/alea_detail.py
diff --git a/lib/pyalps/apptest.py b/bindings/python/pyalps/src/pyalps/apptest.py
similarity index 99%
rename from lib/pyalps/apptest.py
rename to bindings/python/pyalps/src/pyalps/apptest.py
index 02476c784..30389bbc0 100644
--- a/lib/pyalps/apptest.py
+++ b/bindings/python/pyalps/src/pyalps/apptest.py
@@ -485,8 +485,9 @@ def checkProperties( testfile, reffile ):
         del tprop['filename']
         del rprop['filename']
 
-    if cmp(tprop, rprop) == 0:    return True
-    else:   return False
+    if tprop.keys() != rprop.keys():
+        return False
+    return all(np.array_equal(tprop[key], rprop[key]) for key in tprop)
 
 def compareTest( testinputfile, outputs, tmpdir, tstart, compMethod='auto' ):
 
@@ -798,4 +799,3 @@ def createTest( script, inputs=None, outputs=None, prefix=None, refdir='./ref' )
 
     f.close()
     os.chmod(scriptname_prefixed, 0o755)
-
diff --git a/lib/pyalps/cxx.py b/bindings/python/pyalps/src/pyalps/cxx.py
similarity index 80%
rename from lib/pyalps/cxx.py
rename to bindings/python/pyalps/src/pyalps/cxx.py
index dff6bdfc2..6265f2133 100644
--- a/lib/pyalps/cxx.py
+++ b/bindings/python/pyalps/src/pyalps/cxx.py
@@ -33,28 +33,30 @@
 ## while testing (absolute modules, available via PYTHONPATH)
 
 try:
-    from . import pyalea_c
-    from . import pymcdata_c
-    from . import pyngsapi_c
-    from . import pyngsbase_c
-    from . import pyngshdf5_c
-    from . import pyngsobservable_c
-    from . import pyngsobservables_c
-    from . import pyngsparams_c
-    from . import pyngsrandom01_c
-    from . import pyngsresult_c
-    from . import pyngsresults_c
-    from . import pytools_c
+    from ._ext import pyalea_c
+    from ._ext import pymcdata_c
+    from ._ext import pyngsbase_c
+    from ._ext import pyngsapi_c
+    from ._ext import pyngshdf5_c
+    from ._ext import pyngsobservable_c
+    from ._ext import pyngsobservables_c
+    from ._ext import pyngsparams_c
+    from ._ext import pyngsrandom01_c
+    from ._ext import pyngsaccumulator_c
+    from ._ext import pyngsresult_c
+    from ._ext import pyngsresults_c
+    from ._ext import pytools_c
 except ImportError:
     import pyalea_c
     import pymcdata_c
-    import pyngsapi_c
     import pyngsbase_c
+    import pyngsapi_c
     import pyngshdf5_c
     import pyngsobservable_c
     import pyngsobservables_c
     import pyngsparams_c
     import pyngsrandom01_c
+    import pyngsaccumulator_c
     import pyngsresult_c
     import pyngsresults_c
     import pytools_c
diff --git a/lib/pyalps/dataset.py b/bindings/python/pyalps/src/pyalps/dataset.py
similarity index 100%
rename from lib/pyalps/dataset.py
rename to bindings/python/pyalps/src/pyalps/dataset.py
diff --git a/lib/pyalps/dict_intersect.py b/bindings/python/pyalps/src/pyalps/dict_intersect.py
similarity index 89%
rename from lib/pyalps/dict_intersect.py
rename to bindings/python/pyalps/src/pyalps/dict_intersect.py
index 3a683e808..c80df81d6 100644
--- a/lib/pyalps/dict_intersect.py
+++ b/bindings/python/pyalps/src/pyalps/dict_intersect.py
@@ -28,6 +28,12 @@
 
 import numpy as np
 
+def _values_equal(left, right):
+    try:
+        return bool(np.all(left == right))
+    except (TypeError, ValueError):
+        return False
+
 def dict_intersect(dicts):
     """ computes the intersection of a list of dicts
     
@@ -42,12 +48,8 @@ def dict_intersect(dicts):
         take = True
         val0 = dicts[0][key]
         for idict in dicts:
-            try:
-                if val0 != idict[key]:
-                    take = False
-            except:
-                if np.all(val0 != idict[key]):
-                    take = False
+            if not _values_equal(val0, idict[key]):
+                take = False
         if take:
             ret[key] = dicts[0][key]
     return ret
@@ -62,7 +64,7 @@ def dict_difference(dicts):
         take = True
         val0 = dicts[0][key]
         for idict in dicts:
-            if val0 != idict[key]:
+            if not _values_equal(val0, idict[key]):
                 take = False
         if not take:
             ret.append(key)
diff --git a/lib/pyalps/dwa.py b/bindings/python/pyalps/src/pyalps/dwa.py
similarity index 95%
rename from lib/pyalps/dwa.py
rename to bindings/python/pyalps/src/pyalps/dwa.py
index 2745b6440..44f974696 100644
--- a/lib/pyalps/dwa.py
+++ b/bindings/python/pyalps/src/pyalps/dwa.py
@@ -31,14 +31,8 @@
 import os;
 from . import math;
 import numpy;
-import scipy;
-import matplotlib;
-import matplotlib.pyplot;
 import pyalps;
-try:
-    from .dwa_c import worldlines, bandstructure
-except ImportError:
-    from dwa_c import worldlines, bandstructure
+from ._ext.dwa_c import worldlines, bandstructure
 from functools import reduce
 
 
@@ -106,8 +100,8 @@ def thermalized(h5_outfile, observables, tolerance=0.01, simplified=False, inclu
     timeseries = pyalps.hdf5.archive(h5_outfile, 'r')["/simulation/results/" + observable]['timeseries']['data'];
     mean = timeseries.mean();
 
-    index = scipy.linspace(0, timeseries.size-1, timeseries.size);
-    timeseries = scipy.polyval(scipy.polyfit(index, timeseries, 1), index);  # timeseries get fitted
+    index = numpy.linspace(0, timeseries.size-1, timeseries.size);
+    timeseries = numpy.polyval(numpy.polyfit(index, timeseries, 1), index);  # timeseries get fitted
 
     percentage_increment = (timeseries[-1] - timeseries[0])/mean;
     result = abs(percentage_increment) < tolerance;
@@ -206,6 +200,8 @@ def extract_worldlines(infile, outfile=None):
     return wl;
 
 def show_worldlines(wl=None, reshape=None, at=None, scatter_plot=False, Nmax=20, linewidth=2, linespace=0.1):
+  import matplotlib.pyplot as plt
+
   if wl == None:
     return;
 
@@ -231,15 +227,15 @@ def show_worldlines(wl=None, reshape=None, at=None, scatter_plot=False, Nmax=20,
       wl_coordinates.append([idx,wl_time[idx][idx2]]);
   [wl_coordinates_site, wl_coordinates_time] = numpy.array(wl_coordinates).transpose();
 
-  matplotlib.pyplot.figure(frameon=False);
-  matplotlib.pyplot.xticks(range(wl_idx.size), wl_idx);
-  matplotlib.pyplot.yticks([0,1]);
-  matplotlib.pyplot.xlim(-0.5, wl_idx.size-0.5);
-  matplotlib.pyplot.ylim(-0.05,1.05);
+  plt.figure(frameon=False);
+  plt.xticks(range(wl_idx.size), wl_idx);
+  plt.yticks([0,1]);
+  plt.xlim(-0.5, wl_idx.size-0.5);
+  plt.ylim(-0.05,1.05);
 
   if scatter_plot:
-    matplotlib.pyplot.scatter(wl_coordinates_site, wl_coordinates_time);
-    matplotlib.pyplot.show();
+    plt.scatter(wl_coordinates_site, wl_coordinates_time);
+    plt.show();
     return;
 
   wl_state_segments = [];
@@ -275,23 +271,23 @@ def show_worldlines(wl=None, reshape=None, at=None, scatter_plot=False, Nmax=20,
 
   for wl_n_state_segment in wl_n_state_segments[0]:
     [segment_site, segment_time] = numpy.array(wl_n_state_segment).transpose();  
-    matplotlib.pyplot.plot(segment_site, segment_time, '--k', linewidth=linewidth);
+    plt.plot(segment_site, segment_time, '--k', linewidth=linewidth);
   
   for wl_n_state_segment in wl_n_state_segments[1]:
     [segment_site, segment_time] = numpy.array(wl_n_state_segment).transpose();
-    matplotlib.pyplot.plot(segment_site, segment_time, '-k', linewidth=linewidth);
+    plt.plot(segment_site, segment_time, '-k', linewidth=linewidth);
   
   for n in range(2,Nmax+1):
     for wl_n_state_segment in wl_n_state_segments[n]:
       [segment_site, segment_time] = numpy.array(wl_n_state_segment).transpose();
       for m in range(n):
-        matplotlib.pyplot.plot(segment_site - (m - (n-1)/2.)*linespace, segment_time, '-k', linewidth=2);
+        plt.plot(segment_site - (m - (n-1)/2.)*linespace, segment_time, '-k', linewidth=2);
   
   for wl_vertex_segment in wl_vertex_segments:
     [segment_site, segment_time] = numpy.array(wl_vertex_segment).transpose();
-    matplotlib.pyplot.plot(segment_site, segment_time, '-k', linewidth=linewidth);
+    plt.plot(segment_site, segment_time, '-k', linewidth=linewidth);
   
-  matplotlib.pyplot.show();
+  plt.show();
   return;  
 
 def recursiveRun(cmd, cmd_lang='command_line', follow_up_script=None, end_script=None, n=None, break_if=None, break_elseif=None, write_status=None, loc=None, loc0=None, batch_submit=False, batch_cmd_prefix=None, batch_run_directory=None, batch_run_script='run.script', batch_next_run_script=None, batch_run_now=False, batch_noRun=False):
@@ -512,5 +508,3 @@ def summaryReport(h5_outfile):
 
   
 
-
-
diff --git a/lib/pyalps/fit_wrapper.py b/bindings/python/pyalps/src/pyalps/fit_wrapper.py
similarity index 100%
rename from lib/pyalps/fit_wrapper.py
rename to bindings/python/pyalps/src/pyalps/fit_wrapper.py
diff --git a/lib/pyalps/floatwitherror.py b/bindings/python/pyalps/src/pyalps/floatwitherror.py
similarity index 100%
rename from lib/pyalps/floatwitherror.py
rename to bindings/python/pyalps/src/pyalps/floatwitherror.py
diff --git a/lib/pyalps/hdf5.py b/bindings/python/pyalps/src/pyalps/hdf5.py
similarity index 100%
rename from lib/pyalps/hdf5.py
rename to bindings/python/pyalps/src/pyalps/hdf5.py
diff --git a/lib/pyalps/hlist.py b/bindings/python/pyalps/src/pyalps/hlist.py
similarity index 100%
rename from lib/pyalps/hlist.py
rename to bindings/python/pyalps/src/pyalps/hlist.py
diff --git a/lib/pyalps/lattice.py b/bindings/python/pyalps/src/pyalps/lattice.py
similarity index 100%
rename from lib/pyalps/lattice.py
rename to bindings/python/pyalps/src/pyalps/lattice.py
diff --git a/lib/pyalps/load.py b/bindings/python/pyalps/src/pyalps/load.py
similarity index 100%
rename from lib/pyalps/load.py
rename to bindings/python/pyalps/src/pyalps/load.py
diff --git a/lib/pyalps/math.py b/bindings/python/pyalps/src/pyalps/math.py
similarity index 100%
rename from lib/pyalps/math.py
rename to bindings/python/pyalps/src/pyalps/math.py
diff --git a/lib/pyalps/maxent.py b/bindings/python/pyalps/src/pyalps/maxent.py
similarity index 95%
rename from lib/pyalps/maxent.py
rename to bindings/python/pyalps/src/pyalps/maxent.py
index 8bf79ad14..2a7123634 100644
--- a/lib/pyalps/maxent.py
+++ b/bindings/python/pyalps/src/pyalps/maxent.py
@@ -26,8 +26,4 @@
  #                                                                                 #
  # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
 
-try:
-    from .maxent_c import *
-except ImportError:
-    from maxent_c import *
-    
\ No newline at end of file
+from ._ext.maxent_c import *
diff --git a/lib/pyalps/mpi.py b/bindings/python/pyalps/src/pyalps/mpi.py
similarity index 100%
rename from lib/pyalps/mpi.py
rename to bindings/python/pyalps/src/pyalps/mpi.py
diff --git a/lib/pyalps/mpl_setup_macosx.py b/bindings/python/pyalps/src/pyalps/mpl_setup_macosx.py
similarity index 100%
rename from lib/pyalps/mpl_setup_macosx.py
rename to bindings/python/pyalps/src/pyalps/mpl_setup_macosx.py
diff --git a/lib/pyalps/mpl_setup_qt.py b/bindings/python/pyalps/src/pyalps/mpl_setup_qt.py
similarity index 100%
rename from lib/pyalps/mpl_setup_qt.py
rename to bindings/python/pyalps/src/pyalps/mpl_setup_qt.py
diff --git a/lib/pyalps/mpl_setup_tk.py b/bindings/python/pyalps/src/pyalps/mpl_setup_tk.py
similarity index 100%
rename from lib/pyalps/mpl_setup_tk.py
rename to bindings/python/pyalps/src/pyalps/mpl_setup_tk.py
diff --git a/lib/pyalps/natural_sort.py b/bindings/python/pyalps/src/pyalps/natural_sort.py
similarity index 100%
rename from lib/pyalps/natural_sort.py
rename to bindings/python/pyalps/src/pyalps/natural_sort.py
diff --git a/lib/pyalps/ngs.py b/bindings/python/pyalps/src/pyalps/ngs.py
similarity index 79%
rename from lib/pyalps/ngs.py
rename to bindings/python/pyalps/src/pyalps/ngs.py
index 9a39864e7..f4dcb1815 100644
--- a/lib/pyalps/ngs.py
+++ b/bindings/python/pyalps/src/pyalps/ngs.py
@@ -27,22 +27,14 @@
  #                                                                                 #
  # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
 
-import sys
-
-if sys.version_info[:2] >= (3, 8):
-    from collections.abc import MutableMapping
-else:
-    from collections import MutableMapping
-import types
-
+from collections.abc import MutableMapping
 from .cxx.pyngsparams_c import params
-params.__bases__ = (MutableMapping, ) + params.__bases__
 
 from .cxx.pyngsobservable_c import observable
-class ObservableOperators:
-    def __lshift__(self, other):
-        self.append(other)
-observable.__bases__ = (ObservableOperators, ) + observable.__bases__
+def _observable_lshift(self, other):
+    self.append(other)
+    return self
+observable.__lshift__ = _observable_lshift
 
 class RealObservable:
     def __init__(self, name, binnum = 0):
@@ -59,7 +51,6 @@ def addToObservables(self, observables): #rename this with new ALEA
         observables.createRealVectorObservable(self.name, self.binnum)
 
 from .cxx.pyngsobservables_c import observables
-observables.__bases__ = (MutableMapping, ) + observables.__bases__
 
 from .cxx.pyngsobservable_c import createRealObservable #remove this with new ALEA!
 from .cxx.pyngsobservable_c import createRealVectorObservable #remove this with new ALEA!
@@ -68,7 +59,17 @@ def addToObservables(self, observables): #rename this with new ALEA
 from .cxx.pyngsresult_c import observable2result #remove this with new ALEA!
 
 from .cxx.pyngsresults_c import results
-results.__bases__ = (MutableMapping, ) + results.__bases__
+
+# Boost.Python allowed mutating extension-type base classes after creation.
+# nanobind extension types use a different allocator/deallocator layout, so
+# register them as virtual MutableMapping implementations and copy the mixin
+# methods onto the concrete classes instead.
+for _mapping_type in (params, observables, results):
+    MutableMapping.register(_mapping_type)
+    for _method in ("keys", "values", "items", "get", "pop", "popitem",
+                    "clear", "update", "setdefault", "__eq__", "__ne__"):
+        if not hasattr(_mapping_type, _method):
+            setattr(_mapping_type, _method, getattr(MutableMapping, _method))
 
 from .cxx.pyngsbase_c import mcbase
 
diff --git a/lib/pyalps/plot.py b/bindings/python/pyalps/src/pyalps/plot.py
similarity index 100%
rename from lib/pyalps/plot.py
rename to bindings/python/pyalps/src/pyalps/plot.py
diff --git a/lib/pyalps/plot_core.py b/bindings/python/pyalps/src/pyalps/plot_core.py
similarity index 100%
rename from lib/pyalps/plot_core.py
rename to bindings/python/pyalps/src/pyalps/plot_core.py
diff --git a/lib/pyalps/pyalps_config.py b/bindings/python/pyalps/src/pyalps/pyalps_config.py
similarity index 100%
rename from lib/pyalps/pyalps_config.py
rename to bindings/python/pyalps/src/pyalps/pyalps_config.py
diff --git a/lib/pyalps/pyalps_config.py.in b/bindings/python/pyalps/src/pyalps/pyalps_config.py.in
similarity index 51%
rename from lib/pyalps/pyalps_config.py.in
rename to bindings/python/pyalps/src/pyalps/pyalps_config.py.in
index 599df9a7f..150ef0f5e 100644
--- a/lib/pyalps/pyalps_config.py.in
+++ b/bindings/python/pyalps/src/pyalps/pyalps_config.py.in
@@ -1,2 +1,2 @@
 ALPS_XML_INSTALL_DIR="@CMAKE_INSTALL_PREFIX@/lib/xml"
-ALPS_BIN_INSTALL_DIR="@CMAKE_INSTALL_PREFIX@/bin"
\ No newline at end of file
+ALPS_BIN_INSTALL_DIR="@CMAKE_INSTALL_PREFIX@/bin"
diff --git a/lib/pyalps/pytools.py b/bindings/python/pyalps/src/pyalps/pytools.py
similarity index 100%
rename from lib/pyalps/pytools.py
rename to bindings/python/pyalps/src/pyalps/pytools.py
diff --git a/lib/pyalps/tools.py b/bindings/python/pyalps/src/pyalps/tools.py
similarity index 98%
rename from lib/pyalps/tools.py
rename to bindings/python/pyalps/src/pyalps/tools.py
index ee60690d9..cac8d7a1f 100644
--- a/lib/pyalps/tools.py
+++ b/bindings/python/pyalps/src/pyalps/tools.py
@@ -36,6 +36,7 @@
 import sys
 import glob
 from . import math
+import numpy as np
 import scipy.stats
 import copy
 
@@ -231,7 +232,7 @@ def evaluateQWL(infiles, appname='qwl_evaluate', DELTA_T=None, T_MIN=None, T_MAX
     cmdline += make_list(infiles)
     res = executeCommand(cmdline)
     if res != 0:
-      raise Excpetion("Execution error in evaluateQWL: " + str(res))
+      raise RuntimeError("Execution error in evaluateQWL: " + str(res))
     datasets = []
     for infile in infiles:
       datasets.append([])
@@ -566,9 +567,9 @@ def checkSteadyState(sets=None, outfile=None, observable=None, confidenceInterva
   else:
     ts  = pyalps.loadTimeSeries(outfile, observable);  ### y
     N   = ts.size;
-    idx = scipy.linspace(1, N, N);                     ### x
+    idx = np.linspace(1, N, N);                       ### x
 
-    beta1 = scipy.polyfit(idx, ts, 1)[0];              ### slope
+    beta1 = np.polyfit(idx, ts, 1)[0];                 ### slope
     
     ts_std    = np.std(ts, ddof=1);                          ### unbiased estimate of standard deviation in y
     beta1_std = math.sqrt((12.*ts_std*ts_std)/(N * (N*N-1)));   ### unbiased estimate of standard deviation in slope
@@ -852,7 +853,7 @@ def stringListToList(inList):
         #find number of bracketed items (they come in pairs)
         numbrackets=dum.count('[')
         if numbrackets==0 :
-            unbracketed=map(float,dum.replace('[','').replace(']','').replace(' ','').split(','))
+            unbracketed=list(map(float,dum.replace('[','').replace(']','').replace(' ','').split(',')))
             for q in unbracketed:
                 outList.append([q])
         elif numbrackets>0:
@@ -862,16 +863,16 @@ def stringListToList(inList):
                 startInd=dum.find('[',count)
                 finishInd=dum.find(']', count)
                 if startInd>count:
-                    unbracketed=map(float,(dum[count:startInd-1].replace('[','').replace(']','')\
-                    .replace(' ','').split(',')))
+                    unbracketed=list(map(float,(dum[count:startInd-1].replace('[','').replace(']','')\
+                    .replace(' ','').split(','))))
                     for q in unbracketed:
                         outList.append([q])
-                outList.append(map(float,dum[startInd:finishInd+1].replace('[','').\
-                replace(']','').replace(' ','').split(',')))
+                outList.append(list(map(float,dum[startInd:finishInd+1].replace('[','').\
+                replace(']','').replace(' ','').split(','))))
                 count=finishInd+2
             if len(dum)-count>0:
-                unbracketed=map(float,dum[count:len(dum)].replace('[','').replace(']','')\
-                .replace(' ','').split(','))
+                unbracketed=list(map(float,dum[count:len(dum)].replace('[','').replace(']','')\
+                .replace(' ','').split(',')))
                 for q in unbracketed:
                     outList.append([q])
     else:
@@ -1065,4 +1066,3 @@ def CycleMarkers (data, foreach,
             q.props['line'] = all[key] + '-'
 
     return data
-
diff --git a/config/ALPSConfig.cmake.in b/config/ALPSConfig.cmake.in
index 4436acd40..c7469387f 100644
--- a/config/ALPSConfig.cmake.in
+++ b/config/ALPSConfig.cmake.in
@@ -64,16 +64,6 @@ set(ALPS_BLAS_LIBRARIES "@BLAS_LIBRARIES@")
 set(ALPS_BLAS_LIBRARY "@BLAS_LIBRARY@")
 set(ALPS_MKL_INCLUDE_DIR "@MKL_INCLUDE_DIR@")
 
-# Python
-set(ALPS_HAVE_PYTHON "@ALPS_HAVE_PYTHON@")
-set(ALPS_PYTHON_INTERPRETER "@PYTHON_INTERPRETER@")
-set(ALPS_PYTHON_INCLUDE_DIRS "@PYTHON_INCLUDE_DIRS@")
-set(ALPS_PYTHON_NUMPY_INCLUDE_DIR "@PYTHON_NUMPY_INCLUDE_DIR@")
-set(ALPS_PYTHON_LIBRARY "@PYTHON_LIBRARY@")
-set(ALPS_PYTHON_SITE_PKG "@PYTHON_SITE_PKG@")
-set(ALPS_PYTHON_EXTRA_LIBS "@PYTHON_EXTRA_LIBS@")
-set(ALPS_PYTHON_LINK_FOR_SHARED "@PYTHON_LINK_FOR_SHARED@")
-
 # FFTW
 set(ALPS_FFTW_LIBRARIES "@FFTW_LIBRARIES@")
 set(ALPS_FFTW_INCLUDE_DIR "@FFTW_INCLUDE_DIR@")
@@ -130,4 +120,3 @@ set(ALPS_EXTRA_LIBRARIES "@ALPS_EXTRA_LIBRARIES@")
 # list of ALPS and dependent libraries
 set(ALPS_LIBRARIES alps ${ALPS_Boost_LIBRARIES} ${ALPS_EXTRA_LIBRARIES} CACHE STRING "List of ALPS and dependent libraries." FORCE)
 set(ALPS_FORTRAN_LIBRARIES alps_fortran CACHE STRING "List of ALPS-Fortran library." FORCE)
-
diff --git a/config/FindBoostForALPS.cmake b/config/FindBoostForALPS.cmake
index 7038bb766..8fbed54b7 100644
--- a/config/FindBoostForALPS.cmake
+++ b/config/FindBoostForALPS.cmake
@@ -67,91 +67,6 @@ if(ALPS_USE_SYSTEM_BOOST)
       "Upgrade the system Boost installation or disable ALPS_USE_SYSTEM_BOOST.")
   endif()
 
-  # Save Boost_LIBRARIES now — a second find_package(Boost) call below
-  # (for the Python component) would overwrite this variable.
-  set(_alps_boost_libraries_saved ${Boost_LIBRARIES})
-
-  # Python component: library naming varies by Boost/distro version.
-  # Try python<major><minor> (e.g. python311), then python3, then python.
-  if(ALPS_HAVE_PYTHON)
-    set(_alps_python_component "")
-    set(_alps_python_library "")
-    foreach(_pycomp "python${PYVER}" "python3" "python")
-      find_package(Boost QUIET COMPONENTS ${_pycomp})
-      if(Boost_${_pycomp}_FOUND)
-        # Capture Boost_LIBRARIES right here: after a single-component
-        # find_package it contains exactly that one library path.
-        set(_alps_python_library ${Boost_LIBRARIES})
-        set(_alps_python_component ${_pycomp})
-        break()
-      endif()
-    endforeach()
-
-    # Restore the full library list (overwritten by the python find_package).
-    if(_alps_python_component)
-      message(STATUS "Found system Boost.Python component: ${_alps_python_component}")
-      set(Boost_LIBRARIES ${_alps_boost_libraries_saved} ${_alps_python_library})
-    else()
-      message(WARNING
-        "System Boost.Python library not found (tried python${PYVER}, python3, python). "
-        "Python bindings will be disabled.")
-      set(Boost_LIBRARIES ${_alps_boost_libraries_saved})
-      set(ALPS_HAVE_PYTHON OFF)
-      set(BUILD_BOOST_PYTHON OFF)
-    endif()
-  endif()
-
-  # Set ALPS_HAVE_BOOST_NUMPY for Boost >= 1.63 (when boost::python::numpy
-  # was introduced).
-  if(Boost_VERSION_STRING VERSION_GREATER_EQUAL "1.63.0")
-    set(ALPS_HAVE_BOOST_NUMPY ON)
-  endif()
-
-  # Scenario 3: system Boost 1.63-1.86 + NumPy >= 2.0.
-  # (evaluated below; set the flag early so the numpy lib search is guarded by it)
-  # boost::python::numpy in these versions uses deprecated NumPy C API
-  # removed in NumPy 2.0.  Fall back to boost::python::numeric::array,
-  # which uses only the stable NumPy C API.
-  if(ALPS_HAVE_BOOST_NUMPY AND ALPS_HAVE_PYTHON)
-    EXEC_PYTHON_SCRIPT("import numpy; print(numpy.__version__)" _alps_numpy_ver)
-    message(STATUS "NumPy version: ${_alps_numpy_ver}")
-    if(_alps_numpy_ver VERSION_GREATER_EQUAL "2.0.0" AND
-       Boost_VERSION_STRING VERSION_LESS "1.87.0")
-      message(WARNING
-        "System Boost ${Boost_VERSION_STRING} does not support NumPy >= 2.0 "
-        "(requires Boost >= 1.87). "
-        "Falling back to boost::python::numeric::array. "
-        "Upgrade system Boost to >= 1.87 to silence this warning.")
-      set(ALPS_HAVE_BOOST_NUMPY OFF)
-    endif()
-  endif()
-
-  # Boost.NumPy library: only link when ALPS_HAVE_BOOST_NUMPY is still ON
-  # after the scenario-3 check above.  Library naming mirrors python: try
-  # numpy<major><minor>, numpy3, numpy.
-  if(ALPS_HAVE_BOOST_NUMPY AND ALPS_HAVE_PYTHON)
-    set(_alps_boost_libs_before_numpy ${Boost_LIBRARIES})
-    set(_alps_numpy_lib_found "")
-    foreach(_npcomp "numpy${PYVER}" "numpy3" "numpy")
-      find_package(Boost QUIET COMPONENTS ${_npcomp})
-      if(Boost_${_npcomp}_FOUND)
-        message(STATUS "Found system Boost.NumPy component: ${_npcomp}")
-        # Boost_LIBRARIES is now just the numpy lib — capture and restore.
-        set(_alps_numpy_library ${Boost_LIBRARIES})
-        set(Boost_LIBRARIES ${_alps_boost_libs_before_numpy} ${_alps_numpy_library})
-        set(_alps_numpy_lib_found TRUE)
-        break()
-      endif()
-    endforeach()
-    if(NOT _alps_numpy_lib_found)
-      message(WARNING
-        "System Boost.NumPy library not found (tried numpy${PYVER}, numpy3, numpy). "
-        "Falling back to boost::python::numeric::array.")
-      set(Boost_LIBRARIES ${_alps_boost_libs_before_numpy})
-      set(ALPS_HAVE_BOOST_NUMPY OFF)
-    endif()
-  endif()
-
   # Align Boost_INCLUDE_DIR (singular) used elsewhere in the build.
   set(Boost_INCLUDE_DIR ${Boost_INCLUDE_DIRS})
 
@@ -165,8 +80,6 @@ if(ALPS_USE_SYSTEM_BOOST)
   message(STATUS "Using system Boost ${Boost_VERSION_STRING}")
   message(STATUS "  includes:           ${Boost_INCLUDE_DIRS}")
   message(STATUS "  libraries:          ${Boost_LIBRARIES}")
-  message(STATUS "  ALPS_HAVE_BOOST_NUMPY: ${ALPS_HAVE_BOOST_NUMPY}")
-
   return()
 
 endif() # ALPS_USE_SYSTEM_BOOST
diff --git a/config/FindBoostSrc.cmake b/config/FindBoostSrc.cmake
index a9029a24b..76977602a 100644
--- a/config/FindBoostSrc.cmake
+++ b/config/FindBoostSrc.cmake
@@ -40,10 +40,6 @@ if (NOT DEFINED BUILD_BOOST_SYSTEM)
   set(BUILD_BOOST_SYSTEM TRUE)
 endif (NOT DEFINED BUILD_BOOST_SYSTEM)
 
-if (NOT DEFINED BUILD_BOOST_PYTHON)
-  set(BUILD_BOOST_PYTHON TRUE)
-endif(NOT DEFINED BUILD_BOOST_PYTHON)
-
 if (NOT DEFINED BUILD_BOOST_THREAD)
   set(BUILD_BOOST_THREAD TRUE)
 endif (NOT DEFINED BUILD_BOOST_THREAD)
@@ -95,11 +91,6 @@ if(Boost_INCLUDE_DIR)
   MATH(EXPR Boost_SUBMINOR_VERSION "${Boost_VERSION} % 100")
 endif(Boost_INCLUDE_DIR)
 
-if(Boost_VERSION AND NOT Boost_VERSION LESS 106300)
-  # Boost Numpy is compiled if we have >= 1.63
-  set(ALPS_HAVE_BOOST_NUMPY ON)
-endif(Boost_VERSION AND NOT Boost_VERSION LESS 106300)
-
 if(Boost_ROOT_DIR)
   message(STATUS "Found Boost Source: ${Boost_ROOT_DIR}")
   message(STATUS "Boost Version: ${Boost_MAJOR_VERSION}_${Boost_MINOR_VERSION}_${Boost_SUBMINOR_VERSION}")
@@ -110,19 +101,6 @@ else(Boost_ROOT_DIR)
   message(FATAL_ERROR "Boost Source not Found")
 endif(Boost_ROOT_DIR)
 
-if(BUILD_BOOST_PYTHON)
-  EXEC_PYTHON_SCRIPT ("import numpy; print(numpy.__version__)" numpy_ver)
-  MESSAGE(STATUS "numpy version ${numpy_ver}" )
-  if(${numpy_ver} VERSION_GREATER_EQUAL "2.0.0" AND "${Boost_MAJOR_VERSION}.${Boost_MINOR_VERSION}.${Boost_SUBMINOR_VERSION}" VERSION_LESS "1.87.0" )
-    message(WARNING
-      "Boost ${Boost_MAJOR_VERSION}.${Boost_MINOR_VERSION}.${Boost_SUBMINOR_VERSION} "
-      "does not support NumPy >= 2.0 (requires Boost >= 1.87). "
-      "Falling back to boost::python::numeric::array. "
-      "Upgrade to Boost >= 1.87 to silence this warning.")
-    set(ALPS_HAVE_BOOST_NUMPY OFF)
-  endif()
-endif()
-
 # Avoid auto link of Boost library
 add_definitions(-DBOOST_ALL_NO_LIB=1)
 if(BUILD_SHARED_LIBS)
diff --git a/config/FindPythonMod.cmake b/config/FindPythonMod.cmake
deleted file mode 100644
index cec93bf8a..000000000
--- a/config/FindPythonMod.cmake
+++ /dev/null
@@ -1,275 +0,0 @@
-#  Copyright Olivier Parcollet and Matthias Troyer 2010.
-#   Permission is hereby granted, free of charge, to any person obtaining
-#   a copy of this software and associated documentation files (the “Software”),
-#   to deal in the Software without restriction, including without limitation
-#   the rights to use, copy, modify, merge, publish, distribute, sublicense,
-#   and/or sell copies of the Software, and to permit persons to whom the
-#   Software is furnished to do so, subject to the following conditions:
-#  
-#   The above copyright notice and this permission notice shall be included
-#   in all copies or substantial portions of the Software.
-#  
-#   THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS
-#   OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-#   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-#   AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-#   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
-#   FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
-#   DEALINGS IN THE SOFTWARE.
-
-#
-#  Python settings :
-#
-#  This module checks that :
-#  - the python interpreter is working and version >= 2.6
-#  - it has modules : distutils, numpy, tables, scipy
-#
-#  This module defines the variables
-#  - PYTHON_INTERPRETER : name of the python interpreter
-#  - PYTHON_INCLUDE_DIRS : include for compilation
-#  - PYTHON_NUMPY_INCLUDE_DIR : include for compilation with numpy
-#  - PYTHON_LIBRARY : link flags
-#  - PYTHON_SITE_PKG : path to the standard packages of the python interpreter
-#  - PYTHON_EXTRA_LIBS :  libraries which must be linked in when embedding
-#  - PYTHON_LINK_FOR_SHARED :  linking flags needed when building a shared lib for external modules
-
-message(STATUS "Search for Python")
-
-if (NOT PYTHON_INTERPRETER)
-  #find_program(PYTHON_INTERPRETER NAMES python3 python PATHS $ENV{PATH})
-  if(ALPS_PYTHON_WHEEL OR ALPS_BUILD_LIBS_ONLY)
-    find_package(Python COMPONENTS Interpreter Development.Module REQUIRED)
-    set(PYTHON_LIBRARY Python::Module)
-  else(ALPS_PYTHON_WHEEL OR ALPS_BUILD_LIBS_ONLY)
-    find_package(Python COMPONENTS Interpreter Development REQUIRED)
-    set(PYTHON_LIBRARY Python::Python)
-  endif(ALPS_PYTHON_WHEEL OR ALPS_BUILD_LIBS_ONLY)
-  set(PYTHON_INTERPRETER ${Python_EXECUTABLE})
-message(STATUS "LIBS: ${Python_LIBRARY} ${Python_LIBRARIES}")
-  if (NOT PYTHON_INTERPRETER)
-    set (PYTHON_FOUND FALSE)
-  else(NOT PYTHON_INTERPRETER)
-    set(PYTHON_FOUND TRUE)
-  endif(NOT PYTHON_INTERPRETER)
-else (NOT PYTHON_INTERPRETER)
-  set(PYTHON_FOUND TRUE)
-endif (NOT PYTHON_INTERPRETER)
-
-set(PYTHON_MINIMAL_VERSION 3.9)
-
-if (WIN32)
-  MESSAGE (STATUS "Looking for PythonLibs")
-  find_package(PythonLibs)
-endif (WIN32)
-
-IF (PYTHON_FOUND)
-
-  MESSAGE (STATUS "Python interpreter ${PYTHON_INTERPRETER}")
-  #
-  # The function EXEC_PYTHON_SCRIPT executes the_script in  python interpreter
-  # and set the variable of output_var_name in the calling scope
-  #
-  FUNCTION ( EXEC_PYTHON_SCRIPT the_script output_var_name)
-    EXECUTE_PROCESS(COMMAND ${PYTHON_INTERPRETER} -c "${the_script}"
-      OUTPUT_VARIABLE res RESULT_VARIABLE returncode OUTPUT_STRIP_TRAILING_WHITESPACE)
-    IF (NOT returncode EQUAL 0)
-      MESSAGE(FATAL_ERROR "The script : ${the_script} \n did not run properly in the Python interpreter. Check your python installation.")
-    ENDIF (NOT returncode EQUAL 0)
-    SET( ${output_var_name} ${res} PARENT_SCOPE)
-  ENDFUNCTION (EXEC_PYTHON_SCRIPT)
-
-  #
-  # Check the interpreter and its version
-  #
-  EXEC_PYTHON_SCRIPT ("import sys, string; print(sys.version.split()[0])" PYTHON_VERSION)
-#  STRING(COMPARE GREATER ${PYTHON_MINIMAL_VERSION} ${PYTHON_VERSION} PYTHON_VERSION_NOT_OK)
-#  IF (PYTHON_VERSION_NOT_OK)
-  IF( ${PYTHON_VERSION}  VERSION_LESS  ${PYTHON_MINIMAL_VERSION} )
-    MESSAGE(WARNING "Python intepreter version is ${PYTHON_VERSION} . It should be >= ${PYTHON_MINIMAL_VERSION}")
-    SET(PYTHON_FOUND FALSE)
-  ENDIF ()
-  EXEC_PYTHON_SCRIPT("import sys; print('{}{}'.format(sys.version_info.major,sys.version_info.minor))" PYVER) # e.g. 27, 38
-ENDIF (PYTHON_FOUND)
-
-IF (PYTHON_FOUND)
-  if(PYTHON_VERSION VERSION_LESS "3.11")
-  EXEC_PYTHON_SCRIPT ("import distutils " nulle) # check that distutils is there...
-  else()
-  EXEC_PYTHON_SCRIPT ("import sysconfig " nulle) # check that distutils is there...
-  endif()
-  EXEC_PYTHON_SCRIPT ("import numpy" nulle) # check that numpy is there...
-  #EXEC_PYTHON_SCRIPT ("import scipy" nulle) # check that scipy is there...
-  #EXEC_PYTHON_SCRIPT ("import tables" nulle) # check that tables is there...
-  MESSAGE(STATUS "Python interpreter ok : version ${PYTHON_VERSION}" )
-  
-  #
-  # Python function to normalize linker flags
-  #
-  # Goal: CMake has two requiriments on the library flags:
-  # 1. the string cannot start with a spaces
-  # 2. if the string starts with a slash, the argument is interpreted as *a single library name* or a list of libraries
-  #    this is broken if the linker flags are, e.g. "/path/to/lib -framework MyFramework -sysroot /"
-  # --> we need to split the string into a list of elements starting with "/" or "-".
-  # TODO: there might be problems if some path contains spaces
-  set(PYFUNC_NORMALIZE_FLAGS "def normalize_flags(flags):\n flags=flags.strip()\n if flags[0]=='-':return flags\n parts=flags.split(' ', 1)\n if len(parts)>0:return parts[0].strip()+';'+normalize_flags(parts[1])\n return parts[0].strip()\n")
-  
-  #
-  # Check for Python include path
-  #
-  if(PYTHON_VERSION VERSION_LESS "3.11")
-  EXEC_PYTHON_SCRIPT ("import distutils ; from distutils.sysconfig import * ; print(distutils.sysconfig.get_python_inc())"  PYTHON_INCLUDE_DIRS )
-  else()
-  EXEC_PYTHON_SCRIPT ("import sysconfig ; print(sysconfig.get_path('include'))"  PYTHON_INCLUDE_DIRS )
-  endif()
-  message(STATUS "PYTHON_INCLUDE_DIRS =  ${PYTHON_INCLUDE_DIRS}" )
-  mark_as_advanced(PYTHON_INCLUDE_DIRS)
-  FIND_PATH(TEST_PYTHON_INCLUDE patchlevel.h PATHS ${PYTHON_INCLUDE_DIRS} NO_DEFAULT_PATH)
-  if (NOT TEST_PYTHON_INCLUDE)
-    message (ERROR "The Python header files have not been found. Please check that you installed the Python headers and not only the interpreter.")
-  endif (NOT TEST_PYTHON_INCLUDE)
-
-  #
-  # include files for numpy
-  #
-  EXEC_PYTHON_SCRIPT ("import numpy;print(numpy.get_include())" PYTHON_NUMPY_INCLUDE_DIR)
-  MESSAGE(STATUS "PYTHON_NUMPY_INCLUDE_DIR = ${PYTHON_NUMPY_INCLUDE_DIR}" )
-  mark_as_advanced(PYTHON_NUMPY_INCLUDE_DIR)
-
-  #
-  # Check for site packages
-  #
-  if(PYTHON_VERSION VERSION_LESS "3.11")
-  EXEC_PYTHON_SCRIPT ("from distutils.sysconfig import * ;print(get_python_lib(0,0))"
-              PYTHON_SITE_PKG)
-  else()
-  EXEC_PYTHON_SCRIPT ("import sysconfig ; print(sysconfig.get_path('purelib'))"
-              PYTHON_SITE_PKG)
-  endif()
-  MESSAGE(STATUS "PYTHON_SITE_PKG = ${PYTHON_SITE_PKG}" )
-  mark_as_advanced(PYTHON_SITE_PKG)
-    if (NOT WIN32)
-    if(NOT PYTHON_LIBRARY)
-      #
-      # Check for Python library path
-      #
-      #EXEC_PYTHON_SCRIPT ("import string; from distutils.sysconfig import * ;print string.join(get_config_vars('VERSION'))"  PYTHON_VERSION_MAJOR_MINOR)
-      if(PYTHON_VERSION VERSION_LESS "3.11")
-        EXEC_PYTHON_SCRIPT ("import string; from distutils.sysconfig import *; print(' '.join(get_config_vars('LIBDIR')))" PYTHON_LIBRARY_BASE_PATH)
-        # this is the static libpython which is not always correct. it is better to give precedence to the shared one.
-        # EXEC_PYTHON_SCRIPT ("from distutils.sysconfig import *; print(get_config_vars('LIBRARY')[0])" PYTHON_LIBRARY_BASE_FILE)
-        EXEC_PYTHON_SCRIPT ("from distutils.sysconfig import *; print('libpython{}'.format(' '.join(get_config_vars('VERSION'))))" PYTHON_LIBRARY_BASE_FILE)
-      else()
-        
-        EXEC_PYTHON_SCRIPT ("import string; from sysconfig import *; print(' '.join(get_config_vars('LIBDIR')))" PYTHON_LIBRARY_BASE_PATH)
-        # this is the static libpython which is not always correct. it is better to give precedence to the shared one.
-        # EXEC_PYTHON_SCRIPT ("from distutils.sysconfig import *; print(get_config_vars('LIBRARY')[0])" PYTHON_LIBRARY_BASE_FILE)
-        EXEC_PYTHON_SCRIPT ("from sysconfig import *; print('libpython{}'.format(' '.join(get_config_vars('VERSION'))))" PYTHON_LIBRARY_BASE_FILE)
-      endif()
-      IF(BUILD_SHARED_LIBS)
-        FIND_FILE(PYTHON_LIBRARY NAMES "${PYTHON_LIBRARY_BASE_FILE}.so" PATHS ${PYTHON_LIBRARY_BASE_PATH})
-        IF(NOT PYTHON_LIBRARY)
-          FIND_FILE(PYTHON_LIBRARY NAMES "${PYTHON_LIBRARY_BASE_FILE}m.so" PATHS ${PYTHON_LIBRARY_BASE_PATH})
-        ENDIF(NOT PYTHON_LIBRARY)
-        IF(NOT PYTHON_LIBRARY)
-          FIND_FILE(PYTHON_LIBRARY NAMES "${PYTHON_LIBRARY_BASE_FILE}.a" PATHS ${PYTHON_LIBRARY_BASE_PATH})
-        ENDIF(NOT PYTHON_LIBRARY)
-        IF(NOT PYTHON_LIBRARY)
-          FIND_FILE(PYTHON_LIBRARY NAMES "${PYTHON_LIBRARY_BASE_FILE}m.a" PATHS ${PYTHON_LIBRARY_BASE_PATH})
-        ENDIF(NOT PYTHON_LIBRARY)
-      ELSE(BUILD_SHARED_LIBS)
-        FIND_FILE(PYTHON_LIBRARY NAMES "${PYTHON_LIBRARY_BASE_FILE}.a" PATHS ${PYTHON_LIBRARY_BASE_PATH})
-      ENDIF(BUILD_SHARED_LIBS)
-      IF(NOT PYTHON_LIBRARY)
-        # On Debian/Ubuntu system, libpython*.so is located in /usr/lib/`gcc -print-multiarch`
-        execute_process(COMMAND gcc -print-multiarch OUTPUT_VARIABLE TRIPLES)
-        STRING(REGEX REPLACE "\n" "" TRIPLES ${TRIPLES})
-        FIND_FILE(PYTHON_LIBRARY NAMES "${PYTHON_LIBRARY_BASE_FILE}.so" PATHS "/usr/lib/${TRIPLES}")
-        IF(NOT PYTHON_LIBRARY)
-          FIND_FILE(PYTHON_LIBRARY NAMES "${PYTHON_LIBRARY_BASE_FILE}.a" PATHS "/usr/lib/${TRIPLES}")
-        ENDIF(NOT PYTHON_LIBRARY)
-      ENDIF(NOT PYTHON_LIBRARY)
-    endif(NOT PYTHON_LIBRARY)
-    MESSAGE(STATUS "PYTHON_LIBRARY = ${PYTHON_LIBRARY}" )
-    mark_as_advanced(PYTHON_LIBRARY)
-
-      #
-      # libraries which must be linked in when embedding
-      #
-      if(NOT DEFINED PYTHON_EXTRA_LIBS)
-      if(PYTHON_VERSION VERSION_LESS "3.11")
-        EXEC_PYTHON_SCRIPT ("${PYFUNC_NORMALIZE_FLAGS}from distutils.sysconfig import * ;print( normalize_flags( str(get_config_var('LOCALMODLIBS')) + ' ' + str(get_config_var('LIBS')) + ' ' + str(get_config_var('LDFLAGS')) ))"
-                  PYTHON_EXTRA_LIBS)
-      else()
-        EXEC_PYTHON_SCRIPT ("${PYFUNC_NORMALIZE_FLAGS}from sysconfig import * ;print( normalize_flags( str(get_config_var('LOCALMODLIBS')) + ' ' + str(get_config_var('LIBS')) + ' ' + str(get_config_var('LDFLAGS')) ))"
-                  PYTHON_EXTRA_LIBS)
-      endif()
-      endif()
-      MESSAGE(STATUS "PYTHON_EXTRA_LIBS =${PYTHON_EXTRA_LIBS}" )
-      mark_as_advanced(PYTHON_EXTRA_LIBS)
-
-      #
-      # linking flags needed when embedding (building a shared lib)
-      # To BE RETESTED
-      #
-      if(PYTHON_VERSION VERSION_LESS "3.11")
-        EXEC_PYTHON_SCRIPT ("from distutils.sysconfig import *;print(get_config_var('LINKFORSHARED'))"
-                  PYTHON_LINK_FOR_SHARED)
-      else()
-        EXEC_PYTHON_SCRIPT ("from sysconfig import *;print(get_config_var('LINKFORSHARED'))"
-                  PYTHON_LINK_FOR_SHARED)
-     endif()
-      MESSAGE(STATUS "PYTHON_LINK_FOR_SHARED =  ${PYTHON_LINK_FOR_SHARED}" )
-      mark_as_advanced(PYTHON_LINK_FOR_SHARED)
-    endif(NOT WIN32)
-
-  # Correction on Mac
-  IF(APPLE)
-      SET (PYTHON_LINK_FOR_SHARED -u _PyMac_Error -framework Python)
-      SET (PYTHON_LINK_MODULE -bundle -undefined dynamic_lookup)
-  ELSE(APPLE)
-      SET (PYTHON_LINK_MODULE -shared)
-  ENDIF(APPLE)
-ENDIF (PYTHON_FOUND)
-
-set (PYTHONLIBS_FOUND ${PYTHON_FOUND})
-
-
-EXEC_PYTHON_SCRIPT("import sys; print('{}.{}'.format(sys.version_info.major,sys.version_info.minor))" PYVER) # e.g. 27, 38
-set(ALPS_PYTHON_LIB_DEST_ROOT lib/python${PYVER}/site-packages CACHE PATH "Module install path")
-
-#
-# This function writes down a script to compile f2py modules
-# indeed, one needs to use the f2py of the correct numpy module.
-#
-FUNCTION( WriteScriptToBuildF2pyModule filename fcompiler_desc modulename module_pyf_name filelist )
-  # Copy all the files
-  EXECUTE_PROCESS(COMMAND cp ${CMAKE_CURRENT_SOURCE_DIR}/${module_pyf_name} ${CMAKE_CURRENT_BINARY_DIR} )
-  FOREACH( f ${filelist})
-    EXECUTE_PROCESS(COMMAND cp ${CMAKE_CURRENT_SOURCE_DIR}/${f} ${CMAKE_CURRENT_BINARY_DIR} )
-  ENDFOREACH(f)
-  # write the script that will build the f2py extension
-  SET(filename ${CMAKE_CURRENT_BINARY_DIR}/${filename} )
-  FILE(WRITE ${filename} "import sys\n")
-  FILE(APPEND ${filename} "from numpy.f2py import main\n")
-  FILE(APPEND ${filename} "sys.argv = [''] +'-c --fcompiler=${fcompiler_desc} -m ${modulename} ${modulename}.pyf ${filelist} -llapack'.split()\n")
-  FILE(APPEND ${filename} "main()\n")
-ENDFUNCTION(WriteScriptToBuildF2pyModule)
-
-FUNCTION(PYTHON_ADD_MODULE _NAME )
-  OPTION(PYTHON_ENABLE_MODULE_${_NAME} "Add module ${_NAME}" TRUE)
-  OPTION(PYTHON_MODULE_${_NAME}_BUILD_SHARED "Add module ${_NAME} shared" ${BUILD_SHARED_LIBS})
-
-  IF(PYTHON_ENABLE_MODULE_${_NAME})
-    IF(PYTHON_MODULE_${_NAME}_BUILD_SHARED)
-      SET(PY_MODULE_TYPE MODULE)
-    ELSE(PYTHON_MODULE_${_NAME}_BUILD_SHARED)
-      SET(PY_MODULE_TYPE STATIC)
-      SET_PROPERTY(GLOBAL  APPEND  PROPERTY  PY_STATIC_MODULES_LIST ${_NAME})
-    ENDIF(PYTHON_MODULE_${_NAME}_BUILD_SHARED)
-
-    SET_PROPERTY(GLOBAL  APPEND  PROPERTY  PY_MODULES_LIST ${_NAME})
-    ADD_LIBRARY(${_NAME} ${PY_MODULE_TYPE} ${ARGN})
-#    TARGET_LINK_LIBRARIES(${_NAME} ${PYTHON_LIBRARIES})
-
-  ENDIF(PYTHON_ENABLE_MODULE_${_NAME})
-ENDFUNCTION(PYTHON_ADD_MODULE)
diff --git a/config/UseALPS.cmake b/config/UseALPS.cmake
index 0b4160df0..1e9c5b356 100644
--- a/config/UseALPS.cmake
+++ b/config/UseALPS.cmake
@@ -57,15 +57,6 @@ if(NOT ALPS_USE_FILE_INCLUDED)
   set(BLAS_LIBRARY ${ALPS_BLAS_LIBRARY})
   set(MKL_INCLUDE_DIR ${ALPS_MKL_INCLUDE_DIR})
 
-  # Python
-  set(PYTHON_INTERPRETER ${ALPS_PYTHON_INTERPRETER})
-  set(PYTHON_INCLUDE_DIRS ${ALPS_PYTHON_INCLUDE_DIRS})
-  set(PYTHON_NUMPY_INCLUDE_DIR ${ALPS_PYTHON_NUMPY_INCLUDE_DIR})
-  set(PYTHON_LIBRARY ${ALPS_PYTHON_LIBRARY})
-  set(PYTHON_SITE_PKG ${ALPS_PYTHON_SITE_PKG})
-  set(PYTHON_EXTRA_LIBS ${ALPS_PYTHON_EXTRA_LIBS})
-  set(PYTHON_LINK_FOR_SHARED ${ALPS_PYTHON_LINK_FOR_SHARED})
-
   # FFTW
   set(FFTW_LIBRARIES ${ALPS_FFTW_LIBRARIES})
   set(FFTW_INCLUDE_DIR ${ALPS_FFTW_INCLUDE_DIR})
diff --git a/lib/pyalps/CMakeLists.txt b/lib/pyalps/CMakeLists.txt
deleted file mode 100644
index 2e4672079..000000000
--- a/lib/pyalps/CMakeLists.txt
+++ /dev/null
@@ -1,157 +0,0 @@
-#  Copyright Matthias Troyer 2009 - 2010.
-#   Permission is hereby granted, free of charge, to any person obtaining
-#   a copy of this software and associated documentation files (the “Software”),
-#   to deal in the Software without restriction, including without limitation
-#   the rights to use, copy, modify, merge, publish, distribute, sublicense,
-#   and/or sell copies of the Software, and to permit persons to whom the
-#   Software is furnished to do so, subject to the following conditions:
-#  
-#   The above copyright notice and this permission notice shall be included
-#   in all copies or substantial portions of the Software.
-#  
-#   THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS
-#   OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-#   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-#   AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-#   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
-#   FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
-#   DEALINGS IN THE SOFTWARE.
-
-#
-# python exports
-#
-
-set(ALPS_SHARED_CPPFLAGS PYALPS_EXPORTS=1)
-set(ALPS_STATIC_CPPFLAGS "")
-
-if (ALPS_HAVE_PYTHON AND NOT ALPS_BUILD_LIBS_ONLY)
-
-  set(OLD_SHARED ${BUILD_SHARED_LIBS})
-  set(BUILD_SHARED_LIBS ON)
-  set(PYALEA_SOURCES ../../src/alps/python/pyalea.cpp )
-  set(PYMCDATA_SOURCES ../../src/alps/python/pymcdata.cpp )
-  set(PYTOOLS_SOURCES ../../src/alps/python/pytools.cpp)
-
-  set(PYALPS_SOURCES pyalea_c pymcdata_c pytools_c pyngsparams_c pyngshdf5_c pyngsbase_c  
-      pyngsobservable_c pyngsobservables_c pyngsresult_c pyngsresults_c pyngsapi_c pyngsrandom01_c
-  )
-
-  if(ALPS_NGS_USE_NEW_ALEA)
-    list(APPEND PYALPS_SOURCES pyngsaccumulator_c)
-  endif(ALPS_NGS_USE_NEW_ALEA)
-
-  set (MAXENT_SOURCES ../../tool/maxent.cpp ../../tool/maxent_helper.cpp
-                      ../../tool/maxent_simulation.cpp ../../tool/maxent_parms.cpp)
-
-  set (CTHYB_SOURCES ../../applications/dmft/qmc/hybridization/hybmain.cpp
-       ../../applications/dmft/qmc/hybridization/hybsim.cpp
-       ../../applications/dmft/qmc/hybridization/hyblocal.cpp
-       ../../applications/dmft/qmc/hybridization/hybint.cpp
-       ../../applications/dmft/qmc/hybridization/hybfun.cpp
-       ../../applications/dmft/qmc/hybridization/hybretintfun.cpp
-       ../../applications/dmft/qmc/hybridization/hybmatrix.cpp
-       ../../applications/dmft/qmc/hybridization/hybmatrix_ft.cpp
-       ../../applications/dmft/qmc/hybridization/hybconfig.cpp
-       ../../applications/dmft/qmc/hybridization/hybupdates.cpp
-       ../../applications/dmft/qmc/hybridization/hybevaluate.cpp
-       ../../applications/dmft/qmc/hybridization/hybmeasurements.cpp)
-  set (CTINT_SOURCES ../../applications/dmft/qmc/interaction_expansion2/main.cpp
-       ../../applications/dmft/qmc/fouriertransform.C
-       ../../applications/dmft/qmc/interaction_expansion2/auxiliary.cpp
-       ../../applications/dmft/qmc/interaction_expansion2/observables.cpp 
-       ../../applications/dmft/qmc/interaction_expansion2/fastupdate.cpp
-       ../../applications/dmft/qmc/interaction_expansion2/selfenergy.cpp 
-       ../../applications/dmft/qmc/interaction_expansion2/solver.cpp 
-       ../../applications/dmft/qmc/interaction_expansion2/io.cpp 
-       ../../applications/dmft/qmc/interaction_expansion2/splines.cpp 
-       ../../applications/dmft/qmc/interaction_expansion2/interaction_expansion.cpp 
-       ../../applications/dmft/qmc/interaction_expansion2/measurements.cpp 
-       ../../applications/dmft/qmc/interaction_expansion2/model.cpp)
-
-  set(PYNGSPARAMS_SOURCES ../../src/alps/ngs/python/params.cpp)
-  set(PYNGSHDF5_SOURCES ../../src/alps/ngs/python/hdf5.cpp)
-  set(PYNGSBASE_SOURCES ../../src/alps/ngs/python/mcbase.cpp)
-  set(PYNGSOBSERVABLE_SOURCES ../../src/alps/ngs/python/observable.cpp)
-  set(PYNGSOBSERVABLES_SOURCES ../../src/alps/ngs/python/observables.cpp)
-  set(PYNGSRESULT_SOURCES ../../src/alps/ngs/python/result.cpp)
-  set(PYNGSRESULTS_SOURCES ../../src/alps/ngs/python/results.cpp)
-  set(PYNGSAPI_SOURCES ../../src/alps/ngs/python/api.cpp)
-  set(PYNGSRANDOM01_SOURCES ../../src/alps/ngs/python/random01.cpp)
-  set(PYNGSACCUMULATOR_SOURCES ../../src/alps/ngs/python/accumulator.cpp)
-
-  if(LAPACK_FOUND AND ALPS_BUILD_APPLICATIONS)
-    set(PYALPS_SOURCES ${PYALPS_SOURCES} maxent_c dwa_c cthyb ctint)
-    python_add_module(maxent_c ${MAXENT_SOURCES})
-    python_add_module(cthyb ${CTHYB_SOURCES})
-    python_add_module(ctint ${CTINT_SOURCES})
-    python_add_module(dwa_c ../../applications/qmc/dwa/python/dwa.cpp)
-    include_directories(../../applications/qmc/dwa)
-    include_directories(../../applications/dmft/qmc)
-    set_target_properties(maxent_c PROPERTIES COMPILE_FLAGS "-DBUILD_PYTHON_MODULE") 
-    set_target_properties(cthyb PROPERTIES COMPILE_FLAGS "-DBUILD_PYTHON_MODULE") 
-    set_target_properties(ctint PROPERTIES COMPILE_FLAGS "-DBUILD_PYTHON_MODULE") 
-  endif(LAPACK_FOUND AND ALPS_BUILD_APPLICATIONS)
-
-  python_add_module(pyalea_c ${PYALEA_SOURCES})
-  python_add_module(pymcdata_c ${PYMCDATA_SOURCES})
-  python_add_module(pytools_c ${PYTOOLS_SOURCES})
-  python_add_module(pyngsparams_c ${PYNGSPARAMS_SOURCES})
-  python_add_module(pyngshdf5_c ${PYNGSHDF5_SOURCES})
-  python_add_module(pyngsbase_c ${PYNGSBASE_SOURCES})
-  python_add_module(pyngsobservable_c ${PYNGSOBSERVABLE_SOURCES})
-  python_add_module(pyngsobservables_c ${PYNGSOBSERVABLES_SOURCES})
-  python_add_module(pyngsresult_c ${PYNGSRESULT_SOURCES})
-  python_add_module(pyngsresults_c ${PYNGSRESULTS_SOURCES})
-  python_add_module(pyngsapi_c ${PYNGSAPI_SOURCES})
-  python_add_module(pyngsrandom01_c ${PYNGSRANDOM01_SOURCES})
-
-  if(ALPS_NGS_USE_NEW_ALEA)
-    python_add_module(pyngsaccumulator_c ${PYNGSACCUMULATOR_SOURCES})
-  endif(ALPS_NGS_USE_NEW_ALEA)
-
-
-
-  FOREACH (name ${PYALPS_SOURCES})
-    if(BUILD_SHARED_LIBS)
-      set_target_properties(${name} PROPERTIES COMPILE_DEFINITIONS "${ALPS_SHARED_CPPFLAGS}")
-      set_target_properties(${name} PROPERTIES PREFIX "")
-      if(WIN32 AND NOT UNIX)
-        set_target_properties(${name} PROPERTIES SUFFIX ".pyd")
-      endif(WIN32 AND NOT UNIX)
-    endif(BUILD_SHARED_LIBS)
-    target_link_libraries(${name} ${LINK_LIBRARIES} ${BLAS_LIBRARY} ${LAPACK_LIBRARY} ${LAPACK_LINKER_FLAGS})
-    if(ALPS_PYTHON_WHEEL)
-      target_link_libraries(${name} alps_python)
-      if(APPLE)
-        set_target_properties(${name} PROPERTIES INSTALL_RPATH "@loader_path/lib" )
-      else(APPLE)
-        set_target_properties(${name} PROPERTIES INSTALL_RPATH "$ORIGIN/lib" )
-      endif(APPLE)
-    else()
-      target_link_libraries(${name} alps)
-      if(APPLE)
-        set_target_properties(${name} PROPERTIES INSTALL_RPATH "@loader_path/../../.." )
-      else(APPLE)
-        set_target_properties(${name} PROPERTIES INSTALL_RPATH "$ORIGIN/../../.." )
-      endif(APPLE)
-    endif()
-  ENDFOREACH(name)
-
-  #######################################################################
-  # install
-  #######################################################################
-  if(NOT ALPS_PYTHON_WHEEL)
-      install(TARGETS ${PYALPS_SOURCES}
-            COMPONENT python
-            RUNTIME DESTINATION ${ALPS_PYTHON_LIB_DEST_ROOT}/pyalps/bin
-            ARCHIVE DESTINATION ${ALPS_PYTHON_LIB_DEST_ROOT}/pyalps
-            LIBRARY DESTINATION ${ALPS_PYTHON_LIB_DEST_ROOT}/pyalps)
-    else(NOT ALPS_PYTHON_WHEEL)
-      install(TARGETS ${PYALPS_SOURCES}
-            COMPONENT python
-            RUNTIME DESTINATION pyalps/bin
-            ARCHIVE DESTINATION pyalps
-            LIBRARY DESTINATION pyalps)
-  endif(NOT ALPS_PYTHON_WHEEL)
-  set(BUILD_SHARED_LIBS ${OLD_SHARED})
-endif (ALPS_HAVE_PYTHON AND NOT ALPS_BUILD_LIBS_ONLY)
diff --git a/pyproject.toml b/pyproject.toml
deleted file mode 100644
index 17faed7e6..000000000
--- a/pyproject.toml
+++ /dev/null
@@ -1,85 +0,0 @@
-[build-system]
-requires = ["scikit-build-core", "numpy", "scipy"]
-build-backend = "scikit_build_core.build"
-
-[tool.scikit-build]
-wheel.packages = ["python/pyalps"]
-build.verbose = true
-#build-dir = "./build_alps"
-#build.tool-args = ["-j8", "-l13"]
-logging.level = "DEBUG"
-
-[tool.scikit-build.cmake.define]
-Boost_SRC_DIR = {env="Boost_SRC_DIR"}
-CMAKE_CXX_FLAGS = "-fPIC -fpermissive -DBOOST_NO_AUTO_PTR -DBOOST_FILESYSTEM_NO_CXX20_ATOMIC_REF -DBOOST_TIMER_ENABLE_DEPRECATED"
-CMAKE_C_FLAGS = "-fPIC"
-ALPS_PYTHON_WHEEL = "ON"
-ALPS_BUILD_FORTRAN = "ON"
-
-[[tool.scikit-build.overrides]]
-if.platform-system = "^darwin"
-cmake.define.CMAKE_CXX_FLAGS = "-fPIC -stdlib=libc++ -fpermissive -DBOOST_NO_AUTO_PTR -DBOOST_FILESYSTEM_NO_CXX20_ATOMIC_REF -DBOOST_TIMER_ENABLE_DEPRECATED"
-cmake.define.Boost_SRC_DIR = {env="Boost_SRC_DIR"}
-cmake.define.ALPS_BUILD_FORTRAN = "ON"
-cmake.define.CMAKE_C_FLAGS = "-fPIC"
-cmake.define.ALPS_PYTHON_WHEEL = "ON"
-
-[project]
-name = "pyalps"
-version = "2.3.4b1"
-authors = [
-  { name="Sergei Iskakov", email="siskakov@umich.edu" },
-  { name="Fei Lin", email="feilin.physics@gmail.com" }
-]
-license = {text = "MIT License"}
-
-dependencies = ["numpy", "scipy"]
-
-description = "Python Applications and Libraries for Physics Simulations"
-readme = "README-py.md"
-requires-python = ">=3.9"
-classifiers = [
-    "Development Status :: 5 - Production/Stable",
-    'Intended Audience :: Science/Research',
-    'Intended Audience :: Developers',
-    'Programming Language :: C++',
-    'Programming Language :: Python',
-    'Programming Language :: Python :: 3',
-    'Programming Language :: Python :: 3.9',
-    'Programming Language :: Python :: 3.10',
-    'Programming Language :: Python :: 3.11',
-    'Programming Language :: Python :: 3.12',
-    'Programming Language :: Python :: 3.13',
-    'Programming Language :: Python :: 3 :: Only',
-    'Programming Language :: Python :: Implementation :: CPython',
-    "Operating System :: POSIX",
-    "Operating System :: Unix",
-    "Operating System :: MacOS",
-]
-
-[project.urls]
-Homepage = "https://alps.comp-phys.com"
-Issues = "https://github.com/ALPSim/ALPS/issues"
-
-
-[project.optional-dependencies]
-tests = [
-  'coverage>=5.0.3',
-  'pytest',
-  'pytest-benchmark[histogram]>=3.2.1',
-]
-
-[tool.cibuildwheel]
-skip = ["*-musllinux*"]
-test-requires = "pytest"
-test-command = "pytest {project}/test"
-manylinux-x86_64-image = "manylinux_2_28"
-
-[tool.cibuildwheel.linux]
-before-all = "dnf install -y epel-release; dnf config-manager --set-enabled powertools; dnf install -y fftw-devel hdf5-devel openblas-devel wget; pipx install patchelf==0.14.5.0 --force"
-#before-build=""
-test-command = "pytest {project}/test/pyalps"
-
-[tool.cibuildwheel.macos]
-before-all = "brew reinstall hdf5 fftw gfortran"
-test-command = "pytest {project}/test/pyalps"
diff --git a/requirements.txt b/requirements.txt
deleted file mode 100644
index 92825131a..000000000
--- a/requirements.txt
+++ /dev/null
@@ -1,2 +0,0 @@
-numpy<2.1
-scipy
\ No newline at end of file
diff --git a/src/alps/CMakeLists.txt b/src/alps/CMakeLists.txt
index 9c495dc88..462986ae3 100644
--- a/src/alps/CMakeLists.txt
+++ b/src/alps/CMakeLists.txt
@@ -22,7 +22,6 @@
 #
 
 set(ALPS_SOURCES "")
-set(ALPS_PYTHON_SOURCES "")
 set(ALPS_SHARED_CPPFLAGS ALPS_EXPORTS=1)
 set(ALPS_STATIC_CPPFLAGS "")
 
@@ -74,10 +73,6 @@ set(ALPS_SOURCES
     set(ALPS_SOURCES ${ALPS_SOURCES} osiris/xdr.c osiris/xdr_array.c osiris/xdr_float.c osiris/xdr_stdio.c)
   endif (NOT ALPS_HAVE_RPC_XDR_H)
   
-  if (ALPS_BUILD_PYTHON)
-    set(ALPS_PYTHON_SOURCES ${ALPS_PYTHON_SOURCES} ngs/lib/get_numpy_type.cpp hdf5/python.cpp python/numpy_array.cpp)
-  endif (ALPS_BUILD_PYTHON)
-
   # OpenMPI ULFM
   if (ALPS_NGS_OPENMPI_ULFM)
     set(ALPS_SOURCES ${ALPS_SOURCES} ngs/lib/ulfm.cpp)
@@ -100,34 +95,12 @@ set(ALPS_SOURCES ${ALPS_SOURCES}
     ngs/lib/clone.cpp ngs/lib/clone_info.cpp ngs/lib/job.cpp ngs/lib/parapack.cpp ngs/lib/worker_factory.cpp
 )
 
-if(ALPS_HAVE_PYTHON)
-  if(ALPS_PYTHON_WHEEL)
-    set(ALPS_PYTHON_SOURCES ${ALPS_PYTHON_SOURCES} ${ALPS_SOURCES})
-    add_library(alps_python ${ALPS_PYTHON_SOURCES})
-    target_compile_definitions(alps_python PRIVATE ALPS_HAVE_PYTHON)
-    get_target_property(XXX alps_python COMPILE_DEFINITIONS)
-    message(STATUS "ALPS_PYTHON: ${XXX}")
-    if(ALPS_HAVE_BOOST_NUMPY)
-      target_compile_definitions(alps_python INTERFACE ALPS_HAVE_BOOST_NUMPY)
-    endif()
-  else()
-    set(ALPS_SOURCES ${ALPS_PYTHON_SOURCES} ${ALPS_SOURCES})
-  endif()
-endif()
-
 add_library(alps ${ALPS_SOURCES})
 
 find_package(HDF5)
 if (Boost_FOUND)
   # link to ${Boost_LIBRARIES} when precompiled Boost libraries found
   set(ALPS_LINK_LIBS ${Boost_LIBRARIES} ${HDF5_LIBRARIES}) # ${SZIP_LIBRARIES})
-  if(PYTHONLIBS_FOUND)
-    if(ALPS_PYTHON_WHEEL AND PYTHON_VERSION VERSION_GREATER_EQUAL "3.13" AND APPLE)
-      set(ALPS_LINK_LIBS ${ALPS_LINK_LIBS} ${PYTHON_LIBRARY}) # ${PYTHON_EXTRA_LIBS})
-    else(ALPS_PYTHON_WHEEL AND PYTHON_VERSION VERSION_GREATER_EQUAL "3.13" AND APPLE)
-      set(ALPS_LINK_LIBS ${ALPS_LINK_LIBS} ${PYTHON_LIBRARY} ${PYTHON_EXTRA_LIBS})
-    endif(ALPS_PYTHON_WHEEL AND PYTHON_VERSION VERSION_GREATER_EQUAL "3.13" AND APPLE)
-  endif(PYTHONLIBS_FOUND)
   if(MPI_FOUND)
     set(ALPS_LINK_LIBS ${ALPS_LINK_LIBS} ${MPI_LIBRARIES})
     if(MPI_EXTRA_LIBRARY)
@@ -138,13 +111,6 @@ if (Boost_FOUND)
 else (Boost_FOUND)
   # "boost" target available when Boost libraries are built from source
   target_link_libraries(alps ${ALPS_BOOST_LIBRARY_NAME} ${HDF5_LIBRARIES}) # ${SZIP_LIBRARIES})
-  if(ALPS_HAVE_PYTHON)
-    if(ALPS_PYTHON_WHEEL)
-      target_link_libraries(alps_python ${ALPS_BOOST_LIBRARY_NAME} ${ALPS_BOOST_PYTHON_LIBRARY_NAME} ${HDF5_LIBRARIES}) # ${SZIP_LIBRARIES})
-    else()
-      target_link_libraries(alps ${ALPS_BOOST_PYTHON_LIBRARY_NAME}) # ${SZIP_LIBRARIES})
-    endif()
-  endif()
 endif (Boost_FOUND)
 
 if(BUILD_SHARED_LIBS)
@@ -156,25 +122,8 @@ if(BUILD_SHARED_LIBS)
   endif(HDF5_DEFINITIONS)
   
   set_target_properties(alps PROPERTIES COMPILE_DEFINITIONS "${ALPS_SHARED_CPPFLAGS}")
-  if(ALPS_HAVE_PYTHON)
-    if(ALPS_PYTHON_WHEEL)
-      set_target_properties(alps_python PROPERTIES COMPILE_DEFINITIONS "${ALPS_SHARED_CPPFLAGS}")
-      target_compile_definitions(alps_python PUBLIC ALPS_HAVE_PYTHON)
-      if(ALPS_HAVE_BOOST_NUMPY)
-        target_compile_definitions(alps_python PUBLIC ALPS_HAVE_BOOST_NUMPY)
-      endif()
-    else(ALPS_PYTHON_WHEEL)
-      target_compile_definitions(alps PUBLIC ALPS_HAVE_PYTHON)
-      if(ALPS_HAVE_BOOST_NUMPY)
-        target_compile_definitions(alps PUBLIC ALPS_HAVE_BOOST_NUMPY)
-      endif()
-    endif(ALPS_PYTHON_WHEEL)
-  endif()
 else(BUILD_SHARED_LIBS)
   set_target_properties(alps PROPERTIES COMPILE_DEFINITIONS "${ALPS_STATIC_CPPFLAGS}")
-  if(ALPS_HAVE_PYTHON)
-    set_target_properties(alps_python PROPERTIES COMPILE_DEFINITIONS "${ALPS_STATIC_CPPFLAGS}")
-  endif()
 endif(BUILD_SHARED_LIBS)
 
 
@@ -185,12 +134,12 @@ if(MSVC)
 endif(MSVC)
 
 # Set soversion for library
-if(NOT WIN32 AND NOT APPLE AND NOT ALPS_PYTHON_WHEEL)
+if(NOT WIN32 AND NOT APPLE)
   set_target_properties(alps
     PROPERTIES
     SOVERSION "${ALPS_VERSION_MAJOR}"
     VERSION "${ALPS_VERSION_MAJOR}.${ALPS_VERSION_MINOR}.${ALPS_VERSION_PATCH}" )
-endif(NOT WIN32 AND NOT APPLE AND NOT ALPS_PYTHON_WHEEL)
+endif(NOT WIN32 AND NOT APPLE)
 
 
 #boost librt linking
@@ -207,7 +156,6 @@ endif()
 #######################################################################
 # install
 #######################################################################
-if(NOT ALPS_PYTHON_WHEEL)
 install(TARGETS alps COMPONENT libraries
         ARCHIVE DESTINATION lib
         LIBRARY DESTINATION lib
@@ -219,13 +167,3 @@ if(ALPS_BUILD_FORTRAN)
           LIBRARY DESTINATION lib
           RUNTIME DESTINATION bin)
 endif(ALPS_BUILD_FORTRAN)
-else ()
-  install(TARGETS alps_python COMPONENT libraries
-          ARCHIVE DESTINATION pyalps/lib
-          LIBRARY DESTINATION pyalps/lib
-          RUNTIME DESTINATION pyalps/bin)
-  install(TARGETS alps COMPONENT libraries
-          ARCHIVE DESTINATION pyalps/lib
-          LIBRARY DESTINATION pyalps/lib
-          RUNTIME DESTINATION pyalps/bin)
-endif()
diff --git a/src/alps/alea/mcanalyze.hpp b/src/alps/alea/mcanalyze.hpp
index e5a36c30a..c8894eac8 100644
--- a/src/alps/alea/mcanalyze.hpp
+++ b/src/alps/alea/mcanalyze.hpp
@@ -341,7 +341,7 @@ typename average_type< typename TimeseriesType::value_type >::type mean(const Ti
   for (typename const_iterator_type<TimeseriesType>::type iter = range_begin(timeseries); iter != range_end(timeseries); ++iter)
     OUT = OUT + *iter;
 
-  return OUT / double(size(timeseries));
+  return OUT / double(alps::size(timeseries));
 }
 
 
@@ -355,7 +355,7 @@ typename average_type< typename TimeseriesType::value_type >::type variance(cons
   using std::pow;
   using alps::numeric::pow;
 
-  if (size(timeseries) < 2) boost::throw_exception(NotEnoughMeasurementsError());
+  if (alps::size(timeseries) < 2) boost::throw_exception(NotEnoughMeasurementsError());
 
   return_type _mean = mean(timeseries);
   return_type OUT;
@@ -366,7 +366,7 @@ typename average_type< typename TimeseriesType::value_type >::type variance(cons
     OUT = OUT + pow(*iter-_mean, 2.);
   }
 
-  return OUT / double(size(timeseries) - 1);
+  return OUT / double(alps::size(timeseries) - 1);
 }
 
 
@@ -380,7 +380,7 @@ mctimeseries< typename average_type< typename TimeseriesType::value_type >::type
   using boost::numeric::operators::operator/;
   using boost::numeric::operators::operator+;
 
-  std::size_t _size = size(timeseries);
+  std::size_t _size = alps::size(timeseries);
   average_type _mean = alps::alea::mean(timeseries);
   average_type _variance = alps::alea::variance(timeseries);
   mctimeseries< average_type > OUT;
@@ -412,7 +412,7 @@ mctimeseries< typename average_type< typename TimeseriesType::value_type >::type
   using boost::numeric::operators::operator/;
   using boost::numeric::operators::operator+;
 
-  std::size_t _size = size(timeseries);
+  std::size_t _size = alps::size(timeseries);
   average_type _mean = mean(timeseries);
   average_type _variance = variance(timeseries);
   mctimeseries< average_type > OUT;
@@ -511,7 +511,7 @@ typename average_type< typename TimeseriesType::value_type >::type error (const
   using alps::numeric::sqrt;
   using boost::numeric::operators::operator/;
 
-  return sqrt( variance(timeseries) / double(size(timeseries)) );
+  return sqrt( variance(timeseries) / double(alps::size(timeseries)) );
 }
 
 
@@ -547,7 +547,7 @@ mctimeseries< typename average_type<typename TimeseriesType::value_type>::type >
   using boost::numeric::operators::operator/;
 
   return_type _running_mean;
-  _running_mean.resize(size(timeseries) );
+  _running_mean.resize(alps::size(timeseries) );
 
   std::partial_sum(range_begin(timeseries), range_end(timeseries), _running_mean.begin(), alps::numeric::plus<average_type, average_type, average_type>() );
 
@@ -565,13 +565,13 @@ mctimeseries< typename average_type<typename TimeseriesType::value_type>::type >
   using boost::numeric::operators::operator/;
 
   mctimeseries<average_type> _reverse_running_mean;
-  _reverse_running_mean.resize(size(timeseries) );
+  _reverse_running_mean.resize(alps::size(timeseries) );
 
   std::partial_sum(static_cast <std::reverse_iterator<typename const_iterator_type<TimeseriesType>::type> > (range_end(timeseries)),
                    static_cast <std::reverse_iterator<typename const_iterator_type<TimeseriesType>::type> > (range_begin(timeseries)),
                    static_cast <std::reverse_iterator<typename iterator_type<TimeseriesType>::type> > (_reverse_running_mean.end() ), alps::numeric::plus<average_type, average_type, average_type>() );
 
-  std::size_t count = size(timeseries);
+  std::size_t count = alps::size(timeseries);
   for (typename iterator_type<TimeseriesType>::type iter = _reverse_running_mean.begin(); iter != _reverse_running_mean.end(); ++iter)
     *iter = *iter / count--;
 
@@ -661,4 +661,3 @@ ALPS_MCANALYZE_IMPLEMENT_OSTREAM(mctimeseries_view)
 #endif
 
 
-
diff --git a/src/alps/config.h.in b/src/alps/config.h.in
index 813cb0a7d..4f1972085 100644
--- a/src/alps/config.h.in
+++ b/src/alps/config.h.in
@@ -136,16 +136,6 @@
 // Define to 1 if you use Xerces C++ XML parser by Apache Software Foundation.
 #cmakedefine ALPS_HAVE_XERCES_PARSER
 
-//
-// Python
-//
-
-// Define to 1 if you have Python on your system.
-//#cmakedefine ALPS_HAVE_PYTHON
-
-// Define to 1 if Boost Numpy (>=1.63) is available
-//#cmakedefine ALPS_HAVE_BOOST_NUMPY
-
 //
 // OpenMP
 //
diff --git a/src/alps/ngs/numeric/vector.hpp b/src/alps/ngs/numeric/vector.hpp
index 3fc6675a8..b55da395a 100644
--- a/src/alps/ngs/numeric/vector.hpp
+++ b/src/alps/ngs/numeric/vector.hpp
@@ -95,36 +95,12 @@ namespace alps {
                 using boost::numeric::operators::operator+;
                 return lhs + rhs;
             }
-            //------------------- operator + with scalar -------------------
-            template<typename T, std::size_t N>
-            std::vector<T> operator + (std::vector<T> arg, T const & scalar) {
-                std::transform(arg.begin(), arg.end(), arg.begin(), boost::lambda::_1 + scalar);
-                return arg;
-            }
-            template<typename T, std::size_t N>
-            std::vector<T> operator + (T const & scalar, std::vector<T> arg) {
-                std::transform(arg.begin(), arg.end(), arg.begin(), scalar + boost::lambda::_1);
-                return arg;
-            }
-
             //------------------- operator - -------------------
             template<typename T, typename U>
             std::vector<T> operator - (std::vector<T> const & lhs, std::vector<U> const & rhs) {
                 using boost::numeric::operators::operator-;
                 return lhs - rhs;
             }
-            //------------------- operator + with scalar -------------------
-            template<typename T, std::size_t N>
-            std::vector<T> operator - (std::vector<T> arg, T const & scalar) {
-                std::transform(arg.begin(), arg.end(), arg.begin(), boost::lambda::_1 + scalar);
-                return arg;
-            }
-            template<typename T, std::size_t N>
-            std::vector<T> operator - (T const & scalar, std::vector<T> arg) {
-                std::transform(arg.begin(), arg.end(), arg.begin(), scalar + boost::lambda::_1);
-                return arg;
-            }
-
             //------------------- operator * vector-vector-------------------
             template<typename T, typename U>
             std::vector<T> operator * (std::vector<T> const & lhs, std::vector<U> const & rhs) {
@@ -141,12 +117,14 @@ namespace alps {
             //------------------- operator + with scalar -------------------
             template<typename T>
             std::vector<T> operator + (T const & scalar, std::vector<T> lhs) {
-                std::transform(lhs.begin(), lhs.end(), lhs.begin(), bind1st(std::plus<T>(), scalar));
+                std::transform(lhs.begin(), lhs.end(), lhs.begin(),
+                               [scalar](T const & value) { return scalar + value; });
                 return lhs;
             }
             template<typename T>
             std::vector<T> operator + (std::vector<T> lhs, T const & scalar) {
-                std::transform(lhs.begin(), lhs.end(), lhs.begin(), bind2nd(std::plus<T>(), scalar));
+                std::transform(lhs.begin(), lhs.end(), lhs.begin(),
+                               [scalar](T const & value) { return value + scalar; });
                 return lhs;
             }
 
diff --git a/src/boost/CMakeLists.txt b/src/boost/CMakeLists.txt
index 30274bb1a..fb8fb09af 100644
--- a/src/boost/CMakeLists.txt
+++ b/src/boost/CMakeLists.txt
@@ -22,10 +22,7 @@
 #
 
 set(BOOST_SOURCES "")
-set(BOOST_PYTHON_SOURCES "")
-set(BOOST_MPI_PYTHON_SOURCES "")
 set(BOOST_LINK_LIBS "")
-set(BOOST_PYTHON_LINK_LIBS "")
 
 # Boost.Date_Time
 if(BUILD_BOOST_DATE_TIME)
@@ -169,40 +166,6 @@ if(BUILD_BOOST_SERIALIZATION)
   add_definitions(-DBOOST_SERIALIZATION_DYN_LINK=1)
 endif(BUILD_BOOST_SERIALIZATION)
 
-# Boost.Python
-if(BUILD_BOOST_PYTHON)
-  if(PYTHON_VERSION GREATER 3 AND Boost_MAJOR_VERSION EQUAL 1 AND Boost_MINOR_VERSION LESS 63)
-    message(WARNING "Python 3 support requires Boost 1.63.0 or newer. Previous versions might build but have sporadic segemtation faults at the end of the execution.")
-  endif()
-  if(ALPS_PYTHON_WHEEL AND PYTHON_VERSION VERSION_GREATER_EQUAL "3.13" AND APPLE)
-    set(BOOST_PYTHON_LINK_LIBS ${BOOST_LINK_LIBS} ${PYTHON_LIBRARY}) # ${PYTHON_EXTRA_LIBS})
-  else(ALPS_PYTHON_WHEEL AND PYTHON_VERSION VERSION_GREATER_EQUAL "3.13" AND APPLE)
-    set(BOOST_PYTHON_LINK_LIBS ${BOOST_LINK_LIBS} ${PYTHON_LIBRARY} ${PYTHON_EXTRA_LIBS})
-  endif(ALPS_PYTHON_WHEEL AND PYTHON_VERSION VERSION_GREATER_EQUAL "3.13" AND APPLE)
-  set(DIRECTORY "${Boost_ROOT_DIR}/libs/python/src")
-  set(SOURCES dict.cpp errors.cpp exec.cpp import.cpp list.cpp long.cpp
-      module.cpp numeric.cpp object_operators.cpp object_protocol.cpp slice.cpp
-      str.cpp tuple.cpp wrapper.cpp converter/arg_to_python_base.cpp
-      converter/builtin_converters.cpp converter/from_python.cpp
-      converter/registry.cpp converter/type_id.cpp object/class.cpp
-      object/enum.cpp object/function.cpp object/function_doc_signature.cpp
-      object/inheritance.cpp object/iterator.cpp object/life_support.cpp
-      object/pickle_support.cpp object/stl_iterator.cpp
-  )
-  if(Boost_MAJOR_VERSION EQUAL 1 AND Boost_MINOR_VERSION GREATER 62 AND ALPS_HAVE_BOOST_NUMPY)
-    set(SOURCES ${SOURCES}
-      numpy/dtype.cpp numpy/matrix.cpp numpy/ndarray.cpp
-      numpy/numpy.cpp numpy/scalars.cpp numpy/ufunc.cpp
-    )
-  endif(Boost_MAJOR_VERSION EQUAL 1 AND Boost_MINOR_VERSION GREATER 62 AND ALPS_HAVE_BOOST_NUMPY)
-  foreach(S ${SOURCES})
-    if(EXISTS ${DIRECTORY}/${S})
-      set(BOOST_PYTHON_SOURCES ${BOOST_PYTHON_SOURCES} ${DIRECTORY}/${S})
-    endif(EXISTS ${DIRECTORY}/${S})
-  endforeach(S)
-  add_definitions(-DBOOST_PYTHON_SOURCE)
-endif(BUILD_BOOST_PYTHON)
-
 # Boost.System
 if(BUILD_BOOST_SYSTEM)
   set(DIRECTORY "${Boost_ROOT_DIR}/libs/system/src")
@@ -261,43 +224,6 @@ if(BUILD_BOOST_THREAD)
 endif(BUILD_BOOST_THREAD)
 
 
-# Boost.MPI Python bindings
-if (BUILD_BOOST_MPI AND BUILD_BOOST_PYTHON AND Boost_ROOT_DIR)
-  set(DIRECTORY "${Boost_ROOT_DIR}/libs/mpi/src/python")
-  set(SOURCES collectives.cpp py_communicator.cpp datatypes.cpp
-      documentation.cpp py_environment.cpp py_nonblocking.cpp py_exception.cpp
-      py_request.cpp skeleton_and_content.cpp status.cpp py_timer.cpp serialize.cpp
-  )
-  foreach(S ${SOURCES})
-    if(EXISTS ${DIRECTORY}/${S})
-      set(BOOST_MPI_PYTHON_SOURCES ${BOOST_MPI_PYTHON_SOURCES} ${DIRECTORY}/${S})
-    endif(EXISTS ${DIRECTORY}/${S})
-  endforeach(S)
-
-  # renmae mpi module to mpi_c
-  set(BOOST_MPI_PYTHON_SOURCES ${BOOST_MPI_PYTHON_SOURCES} mpi/module.cpp)
-
-  if (BOOST_MPI_PYTHON_SOURCES)
-    python_add_module(mpi_c ${BOOST_MPI_PYTHON_SOURCES})
-    if(BUILD_SHARED_LIBS)
-      set_target_properties(mpi_c PROPERTIES COMPILE_DEFINITIONS "${ALPS_SHARED_CPPFLAGS}")
-    if(WIN32 AND NOT UNIX)
-      set_target_properties(mpi_c PROPERTIES SUFFIX ".pyd")
-    endif(WIN32 AND NOT UNIX)
-    endif (BUILD_SHARED_LIBS)
-
-    set_target_properties(mpi_c PROPERTIES PREFIX "")
-    target_link_libraries(mpi_c ${ALPS_BOOST_LIBRARY_NAME} ${ALPS_BOOST_PYTHON_LIBRARY_NAME} ${BOOST_LINK_LIBS} ${BOOST_PYTHON_LINK_LIBS})
-
-    install(TARGETS mpi_c COMPONENT python
-            RUNTIME DESTINATION bin
-            ARCHIVE DESTINATION ${ALPS_PYTHON_LIB_DEST_ROOT}/pyalps
-            LIBRARY DESTINATION ${ALPS_PYTHON_LIB_DEST_ROOT}/pyalps)
-  endif(BOOST_MPI_PYTHON_SOURCES)
-
-endif (BUILD_BOOST_MPI AND BUILD_BOOST_PYTHON AND Boost_ROOT_DIR)
-
-
 #######################################################################
 # install
 #######################################################################
@@ -307,10 +233,6 @@ endif (BUILD_BOOST_MPI AND BUILD_BOOST_PYTHON AND Boost_ROOT_DIR)
 if (NOT Boost_FOUND)
   add_library(${ALPS_BOOST_LIBRARY_NAME} ${BOOST_SOURCES})
   target_link_libraries(${ALPS_BOOST_LIBRARY_NAME} ${BOOST_LINK_LIBS})
-  if(BUILD_BOOST_PYTHON)
-    add_library(${ALPS_BOOST_PYTHON_LIBRARY_NAME} ${BOOST_PYTHON_SOURCES})
-    target_link_libraries(${ALPS_BOOST_PYTHON_LIBRARY_NAME} ${BOOST_PYTHON_LINK_LIBS})
-  endif()
     # Boost.Test
     if(BUILD_BOOST_TEST)
       set(DIRECTORY "${Boost_ROOT_DIR}/libs/test/src")
@@ -390,32 +312,15 @@ if (NOT Boost_FOUND)
       COMMAND ${CMAKE_COMMAND} -E copy ${LIB_NAME} ${PROJECT_BINARY_DIR}/bin)
   endif(MSVC)
 
-  if(NOT ALPS_PYTHON_WHEEL)
   install(TARGETS ${ALPS_BOOST_LIBRARY_NAME} COMPONENT libraries
             RUNTIME DESTINATION bin
             ARCHIVE DESTINATION lib
             LIBRARY DESTINATION lib)
-    if(BUILD_BOOST_PYTHON)
-      install(TARGETS ${ALPS_BOOST_PYTHON_LIBRARY_NAME} COMPONENT libraries
-            RUNTIME DESTINATION bin
-            ARCHIVE DESTINATION lib
-            LIBRARY DESTINATION lib)
-    endif()
-    if (ALPS_INSTALL_BOOST_TEST)
-      install(TARGETS boost_unit_test_framework boost_test_exec_monitor boost_prg_exec_monitor
+  if (ALPS_INSTALL_BOOST_TEST)
+    install(TARGETS boost_unit_test_framework boost_test_exec_monitor boost_prg_exec_monitor
             COMPONENT libraries
             RUNTIME DESTINATION bin
             ARCHIVE DESTINATION lib
             LIBRARY DESTINATION lib)
-    endif(ALPS_INSTALL_BOOST_TEST)
-  else ()
-    install(TARGETS ${ALPS_BOOST_LIBRARY_NAME}  COMPONENT libraries
-            RUNTIME DESTINATION pyalps/bin
-            ARCHIVE DESTINATION pyalps/lib
-            LIBRARY DESTINATION pyalps/lib)
-    install(TARGETS ${ALPS_BOOST_PYTHON_LIBRARY_NAME} COMPONENT libraries
-            RUNTIME DESTINATION pyalps/bin
-            ARCHIVE DESTINATION pyalps/lib
-            LIBRARY DESTINATION pyalps/lib)
-  endif()
+  endif(ALPS_INSTALL_BOOST_TEST)
 endif (NOT Boost_FOUND)
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index bc1c49038..177b33c49 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -16,7 +16,7 @@
 #   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 #   FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 #   DEALINGS IN THE SOFTWARE.
-if(NOT ALPS_PYTHON_WHEEL AND NOT ALPS_BUILD_LIBS_ONLY)
+if(NOT ALPS_BUILD_LIBS_ONLY)
 add_subdirectory(accumulator)
 add_subdirectory(alea)
 add_subdirectory(fixed_capacity)
@@ -30,7 +30,6 @@ add_subdirectory(osiris)
 add_subdirectory(parameter)
 add_subdirectory(parapack)
 add_subdirectory(parser)
-add_subdirectory(pyalps)
 add_subdirectory(random)
 add_subdirectory(utility)
-endif()
\ No newline at end of file
+endif()
diff --git a/test/pyalps/CMakeLists.txt b/test/pyalps/CMakeLists.txt
deleted file mode 100644
index f1c5271e3..000000000
--- a/test/pyalps/CMakeLists.txt
+++ /dev/null
@@ -1,47 +0,0 @@
-#  Copyright Matthias Troyer, Synge Todo and Lukas Gamper 2009 - 2010.
-#   Permission is hereby granted, free of charge, to any person obtaining
-#   a copy of this software and associated documentation files (the “Software”),
-#   to deal in the Software without restriction, including without limitation
-#   the rights to use, copy, modify, merge, publish, distribute, sublicense,
-#   and/or sell copies of the Software, and to permit persons to whom the
-#   Software is furnished to do so, subject to the following conditions:
-#  
-#   The above copyright notice and this permission notice shall be included
-#   in all copies or substantial portions of the Software.
-#  
-#   THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS
-#   OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-#   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-#   AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-#   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
-#   FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
-#   DEALINGS IN THE SOFTWARE.
-
-include_directories(${PROJECT_BINARY_DIR}/src)
-include_directories(${PROJECT_SOURCE_DIR}/src)
-include_directories(${Boost_ROOT_DIR})
-
-#add_executable(loadobs loadobs.cpp)
-#add_dependencies(loadobs alps)
-#target_link_libraries(loadobs alps)
-#add_alps_test(loadobs)
-
-enable_testing()
-if (ALPS_BUILD_PYTHON AND BUILD_SHARED_LIBS)
-#    FOREACH (name pyioarchive pyhdf5io numpylarge pyparams hlist_test mcdata pyhdf5 mcanalyze, accumulators)
-    FOREACH (name pyioarchive pyhdf5io_test pyparams_test hlist_test mcdata_test pyhdf5_test mcanalyze)
-        add_test(python_${name}
-            ${CMAKE_COMMAND}
-                -Dpython_interpreter=${PYTHON_INTERPRETER}
-                -Dcmd=${name}.py
-                -Dinput=${name}
-                -Doutput=${name}
-                -Dpythonpath=${PROJECT_BINARY_DIR}/lib/pyalps:${PROJECT_SOURCE_DIR}/lib
-                -Dsourcedir=${CMAKE_CURRENT_SOURCE_DIR}
-                -Dbinarydir=${CMAKE_CURRENT_BINARY_DIR}
-                -Dcmddir=${CMAKE_CURRENT_SOURCE_DIR}
-                -P ${CMAKE_CURRENT_SOURCE_DIR}/run_python_test.cmake
-        )
-        set_property(TEST python_${name} PROPERTY LABELS pyalps)
-    ENDFOREACH(name)
-ENDIF(ALPS_BUILD_PYTHON AND BUILD_SHARED_LIBS)
diff --git a/test/pyalps/test_binding_surface.py b/test/pyalps/test_binding_surface.py
new file mode 100644
index 000000000..091913c22
--- /dev/null
+++ b/test/pyalps/test_binding_surface.py
@@ -0,0 +1,195 @@
+#!/usr/bin/env python3
+# Copyright (C) 2026 by the ALPS collaboration
+# SPDX-License-Identifier: MIT
+
+"""Lock the public pyalps extension surface after the nanobind migration."""
+
+from __future__ import annotations
+
+import copy
+import importlib
+import os
+import tempfile
+from types import SimpleNamespace
+
+import numpy as np
+
+
+def test_extension_import_surface():
+    import pyalps
+    import pyalps.cxx as cxx
+
+    expected = {
+        "pyalea_c",
+        "pymcdata_c",
+        "pytools_c",
+        "pyngsparams_c",
+        "pyngshdf5_c",
+        "pyngsbase_c",
+        "pyngsobservable_c",
+        "pyngsobservables_c",
+        "pyngsresult_c",
+        "pyngsresults_c",
+        "pyngsapi_c",
+        "pyngsrandom01_c",
+        "pyngsaccumulator_c",
+    }
+    assert pyalps is not None
+    assert expected <= set(vars(cxx))
+
+
+def test_cross_module_parameter_archive_and_rng_roundtrip():
+    from pyalps.cxx import pyngshdf5_c, pyngsparams_c, pyngsrandom01_c
+
+    with tempfile.TemporaryDirectory() as directory:
+        path = os.path.join(directory, "surface.h5")
+        params = pyngsparams_c.params()
+        params["integer"] = 42
+        params["real"] = 3.25
+        params["flag"] = True
+        params["text"] = "nanobind-lock"
+
+        rng = pyngsrandom01_c.random01(91)
+        for _ in range(7):
+            rng()
+
+        archive = pyngshdf5_c.hdf5_archive_impl(path, "w")
+        archive.create_group("/parameters")
+        archive.set_context("/parameters")
+        params.save(archive)
+        archive.set_context("/")
+        rng.save(archive)
+        del archive
+
+        loaded = pyngsparams_c.params()
+        restored_rng = pyngsrandom01_c.random01(0)
+        archive = pyngshdf5_c.hdf5_archive_impl(path, "r")
+        archive.set_context("/parameters")
+        loaded.load(archive)
+        archive.set_context("/")
+        restored_rng.load(archive)
+        del archive
+
+        assert sorted(loaded) == sorted(params)
+        assert int(loaded["integer"]) == 42
+        assert float(loaded["real"]) == 3.25
+        assert bool(loaded["flag"]) is True
+        assert str(loaded["text"]) == "nanobind-lock"
+        assert [rng() for _ in range(5)] == [restored_rng() for _ in range(5)]
+
+
+def test_alea_numpy_and_mcdata_operators():
+    from pyalps.cxx.pyalea_c import MCScalarTimeseries, RealObservable, mean, size
+    from pyalps.cxx.pymcdata_c import MCScalarData
+
+    observable = RealObservable("energy")
+    for sample in (0.9, 1.0, 1.1, 1.0):
+        observable << sample
+    assert observable.count == 4
+    assert abs(observable.mean - 1.0) < 1e-12
+    assert observable.error >= 0
+
+    series = MCScalarTimeseries(np.asarray([1.0, 2.0, 3.0]))
+    assert size(series) == 3
+    assert mean(series) == 2.0
+    np.testing.assert_allclose(series.timeseries(), [1.0, 2.0, 3.0])
+
+    first = MCScalarData(1.0, 0.1)
+    second = MCScalarData(2.0, 0.2)
+    total = first + second
+    assert total.mean == 3.0
+    assert total.error > 0
+    duplicate = copy.deepcopy(total)
+    assert duplicate.mean == total.mean
+    assert duplicate.error == total.error
+
+
+def test_ngs_observable_containers():
+    from pyalps import ngs
+
+    observables = ngs.observables()
+    observables.createRealObservable("magnetization")
+    observables["magnetization"] << 1.5
+    assert "magnetization" in observables
+    assert ngs.observable2result(observables["magnetization"]).count == 1
+
+
+def test_name_encoding_roundtrip():
+    from pyalps.cxx.pytools_c import hdf5_name_decode, hdf5_name_encode
+
+    for value in ("plain", "with space", "slash/inside", "café"):
+        assert hdf5_name_decode(hdf5_name_encode(value)) == value
+
+
+def test_accumulator_surface():
+    from pyalps.cxx.pyngsaccumulator_c import error_accumulator
+
+    accumulator = error_accumulator()
+    for sample in (1.0, 2.0, 3.0):
+        accumulator(sample)
+    result = accumulator.result()
+    assert result.count() == 3
+    assert result.mean() == 2.0
+    assert result.error() >= 0
+
+
+def test_optional_application_extension_surface():
+    for name in ("maxent_c", "dwa_c", "cthyb", "ctint"):
+        module = importlib.import_module("pyalps._ext." + name)
+        assert module.__name__.endswith(name)
+
+
+def test_current_python_numpy_and_scipy_compatibility(monkeypatch):
+    import pyalps
+    import pyalps.dwa as dwa
+
+    assert callable(dwa.thermalized)
+
+    parsed = pyalps.stringListToList("[1,[2,3],4]")
+    assert parsed == [[1.0], [2.0, 3.0], [4.0]]
+
+    shared = pyalps.dict_intersect([
+        {"array": np.array([1, 2]), "scalar": 3},
+        {"array": np.array([1, 2]), "scalar": 3},
+    ])
+    np.testing.assert_array_equal(shared["array"], [1, 2])
+    assert shared["scalar"] == 3
+
+    monkeypatch.setattr(
+        pyalps,
+        "loadTimeSeries",
+        lambda *_args: np.array([1.0, 1.1, 0.9, 1.0]),
+    )
+    steady = pyalps.checkSteadyState(outfile="unused.h5", observable="energy")
+    assert isinstance(steady["value"], (bool, np.bool_))
+
+
+def test_python3_property_comparison(monkeypatch):
+    import pyalps
+    import pyalps.apptest as apptest
+
+    properties = {
+        "test.h5": {"vector": np.array([1, 2]), "value": 3},
+        "reference.h5": {"vector": np.array([1, 2]), "value": 3},
+    }
+
+    class Loader:
+        def GetProperties(self, filenames):
+            return [SimpleNamespace(props=properties[filenames[0]].copy())]
+
+    monkeypatch.setattr(pyalps.load, "Hdf5Loader", Loader)
+    assert apptest.checkProperties("test.h5", "reference.h5")
+
+
+if __name__ == "__main__":
+    for test in (
+        test_extension_import_surface,
+        test_cross_module_parameter_archive_and_rng_roundtrip,
+        test_alea_numpy_and_mcdata_operators,
+        test_ngs_observable_containers,
+        test_name_encoding_roundtrip,
+        test_accumulator_surface,
+        test_optional_application_extension_surface,
+    ):
+        test()
+    print("pyalps binding surface: green")
diff --git a/tool/CMakeLists.txt b/tool/CMakeLists.txt
index 519386439..f90b0f4b4 100644
--- a/tool/CMakeLists.txt
+++ b/tool/CMakeLists.txt
@@ -76,25 +76,14 @@ endif(SQLite_FOUND)
   endif(UNIX AND NOT WIN32)
 
 #
-# lattice-preview and helper program
-#
+  # lattice-preview and helper program
+  #
 
   configure_file(config.py.in ${CMAKE_CURRENT_BINARY_DIR}/config.py)
-  if(WIN32 AND NOT UNIX AND ALPS_BUILD_PYTHON)
-    # in the function add_pi_executable is not present ...
-    option(ALPS_HAS_CMAKE_PI_MACROS "Ignore the PI macros if they are not present" ON)
-    mark_as_advanced(ALPS_HAS_CMAKE_PI_MACROS)
-    if (ALPS_HAS_CMAKE_PI_MACROS)
-      add_pi_executable(lattice-preview preview.py ${CMAKE_CURRENT_BINARY_DIR}/config.py license.py)
-      file(GLOB pi_generated_files ${CMAKE_CURRENT_BINARY_DIR}/lattice-preview/*)
-      install(FILES ${pi_generated_files} DESTINATION bin COMPONENT tools)
-    endif (ALPS_HAS_CMAKE_PI_MACROS)
-  else(WIN32 AND NOT UNIX)
-    configure_file(lattice-preview.in ${CMAKE_CURRENT_BINARY_DIR}/lattice-preview)
-    install(PROGRAMS ${CMAKE_CURRENT_BINARY_DIR}/lattice-preview DESTINATION bin COMPONENT tools)
-    install(FILES ${CMAKE_CURRENT_BINARY_DIR}/config.py preview.py license.py
-      DESTINATION lib/python/alps COMPONENT tools)
-  endif(WIN32 AND NOT UNIX AND ALPS_BUILD_PYTHON)
+  configure_file(lattice-preview.in ${CMAKE_CURRENT_BINARY_DIR}/lattice-preview)
+  install(PROGRAMS ${CMAKE_CURRENT_BINARY_DIR}/lattice-preview DESTINATION bin COMPONENT tools)
+  install(FILES ${CMAKE_CURRENT_BINARY_DIR}/config.py preview.py license.py
+    DESTINATION lib/python/alps COMPONENT tools)
 
   #
   # Analytical continuation with MaxEnt
@@ -114,30 +103,3 @@ endif(SQLite_FOUND)
     target_link_libraries(maxent_linear_grid_numeric alps ${LAPACK_LIBRARY} ${BLAS_LIBRARY})
     add_alps_test(maxent_linear_grid_numeric)
    endif(LAPACK_FOUND)
-
-  #
-  # alpspython script
-  #
-  
-  set(ALPSPYTHON_CONFIGURED FALSE)
-  if(PYTHON_INTERPRETER)
-    if (NOT WIN32)
-      set(PYTHONPATH "${ALPS_PYTHON_LIB_DEST_ROOT}")
-      set(PYTHONBIN "${PYTHON_INTERPRETER}")
-      string(CONFIGURE [[
-        set(PROJECT_SOURCE_DIR "@PROJECT_SOURCE_DIR@")
-        set(PYTHONBIN "@PYTHONBIN@")
-        set(PYTHONPATH "${CMAKE_INSTALL_PREFIX}/@PYTHONPATH@")
-        message(STATUS ": ${PYTHONPATH}")
-        message(STATUS ": ${PYTHONBIN}")
-        configure_file(${PROJECT_SOURCE_DIR}/tool/alpspython.in ${CMAKE_INSTALL_PREFIX}/bin/alpspython )
-        ]] install_script @ONLY)
-      install(CODE ${install_script})
-    else (NOT WIN32)
-      set(PYTHONPATH "%HOMEDRIVE%\\Program Files\\ALPS\\lib;%HOMEDRIVE%\\Program Files (x86)\\ALPS\\lib")
-      set(PYTHONBIN "python")
-      configure_file(alpspython.bat.in ${PROJECT_BINARY_DIR}/tool/alpspython.bat)
-      install(PROGRAMS ${PROJECT_BINARY_DIR}/tool/alpspython.bat DESTINATION bin COMPONENT tools)
-    endif (NOT WIN32)
-    set(ALPSPYTHON_CONFIGURED TRUE)
-  endif(PYTHON_INTERPRETER)
diff --git a/tool/maxent.cpp b/tool/maxent.cpp
index eb43d8d96..b3aaa6af2 100644
--- a/tool/maxent.cpp
+++ b/tool/maxent.cpp
@@ -59,11 +59,11 @@ bool stop_callback(boost::posix_time::ptime const & end_time) {
 
 
 #ifdef BUILD_PYTHON_MODULE
-//compile it as a python module (requires boost::python library)
-using namespace boost::python;
+#include "dict_to_params.hpp"
+namespace nb = nanobind;
 
-void run_it(boost::python::dict parms_){
-    alps::parameters_type<MaxEntSimulation>::type parms(parms_);
+void run_it(nb::dict const & parms_){
+    alps::parameters_type<MaxEntSimulation>::type parms = pyalps::params_from_dict(parms_);
     std::string out_file = boost::lexical_cast<std::string>(parms["BASENAME"]|"results")+std::string(".out.h5");
 
 #else
@@ -104,9 +104,7 @@ void run_it(boost::python::dict parms_){
 }
     
 #ifdef BUILD_PYTHON_MODULE
-    BOOST_PYTHON_MODULE(maxent_c)
-    {
-        def("AnalyticContinuation",run_it);//define python-callable run method
-    };
+    NB_MODULE(maxent_c, m) {
+        m.def("AnalyticContinuation", run_it);
+    }
 #endif
-
diff --git a/tutorials/CMakeLists.txt b/tutorials/CMakeLists.txt
index d88dc6974..bc38361a8 100644
--- a/tutorials/CMakeLists.txt
+++ b/tutorials/CMakeLists.txt
@@ -17,7 +17,7 @@
 #   FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 #   DEALINGS IN THE SOFTWARE.
 
-if(NOT ALPS_PYTHON_WHEEL AND NOT ALPS_BUILD_LIBS_ONLY)
+if(NOT ALPS_BUILD_LIBS_ONLY)
   install(DIRECTORY . DESTINATION tutorials COMPONENT tutorials
           FILES_MATCHING PATTERN "*.py" PATTERN "*.ipynb" PATTERN "*.sh" PATTERN "parm*" PATTERN "*params" PATTERN "*.ip" PATTERN "*.op" PATTERN "*.dat"
           PATTERN "*.parm" PATTERN "*.xml" PATTERN "*input*" PATTERN "*.pvsm"
diff --git a/tutorials/code-07-mcmain-mcbase/CMakeLists.txt b/tutorials/code-07-mcmain-mcbase/CMakeLists.txt
index 5dddd97e8..92dc4a164 100644
--- a/tutorials/code-07-mcmain-mcbase/CMakeLists.txt
+++ b/tutorials/code-07-mcmain-mcbase/CMakeLists.txt
@@ -24,22 +24,3 @@ if (MPI_FOUND)
 	target_link_libraries(mpi_pscan ${ALPS_LIBRARIES})
 
 endif (MPI_FOUND)
-
-if (ALPS_HAVE_PYTHON)
-
-	# rule for generating python export
-	set_property(GLOBAL APPEND PROPERTY PY_MODULES_LIST ising_c)
-	if(BUILD_SHARED_LIBS)
-	    add_library(ising_c MODULE ising.cpp export.cpp)
-		set_target_properties(ising_c PROPERTIES COMPILE_DEFINITIONS "${ALPS_SHARED_CPPFLAGS}")
-		if(WIN32 AND NOT UNIX)
-			set_target_properties(ising_c PROPERTIES SUFFIX ".pyd")
-	    endif(WIN32 AND NOT UNIX)
-	else(BUILD_SHARED_LIBS)
-		set_property(GLOBAL APPEND PROPERTY PY_STATIC_MODULES_LIST ising_c)
-		add_library(ising_c STATIC ising.cpp export.cpp)
-	endif (BUILD_SHARED_LIBS)
-	set_target_properties(ising_c PROPERTIES PREFIX "")
-	target_link_libraries(ising_c ${ALPS_LIBRARIES})
-
-endif (ALPS_HAVE_PYTHON)
diff --git a/tutorials/code-07-mcmain-mcbase/heisenberg/o_n_model/CMakeLists.txt b/tutorials/code-07-mcmain-mcbase/heisenberg/o_n_model/CMakeLists.txt
index 584178bb0..09c3b0f11 100644
--- a/tutorials/code-07-mcmain-mcbase/heisenberg/o_n_model/CMakeLists.txt
+++ b/tutorials/code-07-mcmain-mcbase/heisenberg/o_n_model/CMakeLists.txt
@@ -15,23 +15,3 @@ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${bench_flags}")
 # rule for generating the heisenberg example program
 add_executable(heisenberg heisenberg.cpp)
 target_link_libraries(heisenberg ${ALPS_LIBRARIES})
- 
-
-if (ALPS_HAVE_PYTHON)
-
-	# rule for generating python export
-	set_property(GLOBAL APPEND PROPERTY PY_MODULES_LIST pyndsim)
-	if(BUILD_SHARED_LIBS)
-        add_library(pyndsim MODULE EXCLUDE_FROM_ALL export.cpp)
-		set_target_properties(pyndsim PROPERTIES COMPILE_DEFINITIONS "${ALPS_SHARED_CPPFLAGS}")
-		if(WIN32 AND NOT UNIX)
-			set_target_properties(pyndsim PROPERTIES SUFFIX ".pyd")
-	    endif(WIN32 AND NOT UNIX)
-	else(BUILD_SHARED_LIBS)
-		set_property(GLOBAL APPEND PROPERTY PY_STATIC_MODULES_LIST pyndsim)
-        add_library(pyndsim STATIC EXCLUDE_FROM_ALL export.cpp)
-	endif (BUILD_SHARED_LIBS)
-	set_target_properties(pyndsim PROPERTIES PREFIX "")
-	target_link_libraries(pyndsim ${ALPS_LIBRARIES})
-
-endif (ALPS_HAVE_PYTHON)