updated build system and deprecated modules

.gitignore

@@ -8,3 +8,4 @@ polychrom/_polymer_math.cpp
 *.dat
 .idea
 dist
+polychrom/*.so

docs/conf.py

@@ -63,7 +63,6 @@ def setup(app):
     "simtk.unit",
     "simtk.unit.nanometer",
     "simtk.openmm",
-    "joblib",
     "scipy.interpolate.fitpack2",
 ]
 for mod_name in MOCK_MODULES:

docs/index.rst

@@ -17,7 +17,7 @@ Polychrom requires OpenMM, which can be installed through conda: ``conda install
 
 CUDA is the fastest GPU-assisted backend to OpenMM. You would need to have the required version of CUDA, or install OpenMM compiled for your version of CUDA. 
 
-Other dependencies are simple, and are listed in requirements.txt. All but joblib are installable from either conda/pip, and joblib installs well with pip.
+Other dependencies are simple, and are listed in requirements.txt. All are installable through pip or conda.
 
 Installation errors and possible fixes
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

polychrom/__init__.py

@@ -1 +1,19 @@
 __version__ = "0.1.1"
+
+# Check for OpenMM installation
+try:
+    import openmm
+except ImportError:
+    import warnings
+    warnings.warn(
+        "\n"
+        "OpenMM is not installed. Polychrom requires OpenMM for molecular dynamics simulations.\n"
+        "\n"
+        "Please install OpenMM with the appropriate backend for your system:\n"
+        "  - For CUDA: pip install openmm[cuda13]\n"
+        "  - For CPU:  pip install openmm\n"
+        "\n"
+        "Visit https://openmm.org for more information.",
+        ImportWarning,
+        stacklevel=2
+    )

polychrom/__polymer_math.cpp

@@ -6,7 +6,6 @@
 #include <cmath>
 #include <vector>
 #include <ctime>
-#include <omp.h>
 
 using namespace std;
 

polychrom/cli/show

@@ -5,7 +5,6 @@ import os
 import sys
 import tempfile
 
-import joblib
 import numpy as np
 
 usage = """
@@ -43,7 +42,7 @@ start, end, step will basically select data[start:end:step]
 
 
 if len(sys.argv) < 2:
-    print(useage)
+    print(usage)
     exit()
 
 

polychrom/cli/xyz

@@ -6,7 +6,6 @@ import sys
 import tempfile
 import textwrap
 
-import joblib
 import numpy as np
 
 usage = """

polychrom/legacy/polymerutils.py

@@ -0,0 +1 @@
+from polychrom.polymerutils import *  # noqa: F403

polychrom/polymer_analyses.py

@@ -11,7 +11,7 @@
 ------------------------------
 
 The main function calculating contacts is: :py:func:`polychrom.polymer_analyses.calculate_contacts`
-Right now it is a simple wrapper around scipy.cKDTree.
+Right now it is a simple wrapper around scipy.KDTree.
 
 Another function :py:func:`polychrom.polymer_analyses.smart_contacts` was added recently to help build contact maps
 with a large contact radius. It randomly sub-samples the monomers; by default selecting N/cutoff monomers. It then
@@ -37,7 +37,7 @@
 
 import numpy as np
 import pandas as pd
-from scipy.spatial import cKDTree
+from scipy.spatial import KDTree
 from scipy.ndimage import gaussian_filter1d
 
 try:
@@ -66,7 +66,7 @@ def calculate_contacts(data, cutoff=1.7):
     if np.isnan(data).any():
         raise RuntimeError("Data contains NANs")
 
-    tree = cKDTree(data)
+    tree = KDTree(data)
     pairs = tree.query_pairs(cutoff, output_type="ndarray")
     return pairs
 
@@ -573,7 +573,7 @@ def calculate_cistrans(data, chains, chain_id=0, cutoff=5, pbc_box=False, box_si
     chain_end = chains[chain_id][1]
 
     # all contact pairs available in the scaled data
-    tree = cKDTree(data_scaled, boxsize=box_size)
+    tree = KDTree(data_scaled, boxsize=box_size)
     pairs = tree.query_pairs(cutoff, output_type="ndarray")
 
     # total number of contacts of the marked chain:
@@ -582,7 +582,7 @@ def calculate_cistrans(data, chains, chain_id=0, cutoff=5, pbc_box=False, box_si
     all_signal = len(pairs[pairs < chain_end]) - len(pairs[pairs < chain_start])
 
     # contact pairs of the marked chain with itself
-    tree = cKDTree(data[chain_start:chain_end], boxsize=None)
+    tree = KDTree(data[chain_start:chain_end], boxsize=None)
     pairs = tree.query_pairs(cutoff, output_type="ndarray")
 
     # doubled number of contacts of the marked chain with itself (i.e. cis signal)
@@ -593,3 +593,12 @@ def calculate_cistrans(data, chains, chain_id=0, cutoff=5, pbc_box=False, box_si
     trans_signal = all_signal - cis_signal
 
     return cis_signal, trans_signal
+
+
+def rotation_matrix(rotate):
+    """Calculates rotation matrix based on three rotation angles"""
+    tx, ty, tz = rotate
+    Rx = np.array([[1, 0, 0], [0, np.cos(tx), -np.sin(tx)], [0, np.sin(tx), np.cos(tx)]])
+    Ry = np.array([[np.cos(ty), 0, -np.sin(ty)], [0, 1, 0], [np.sin(ty), 0, np.cos(ty)]])
+    Rz = np.array([[np.cos(tz), -np.sin(tz), 0], [np.sin(tz), np.cos(tz), 0], [0, 0, 1]])
+    return np.dot(Rx, np.dot(Ry, Rz))

polychrom/polymerutils.py

@@ -24,33 +24,23 @@
         xyz = data["pos"]
 """
 
-
 from __future__ import absolute_import, division, print_function, unicode_literals
 
 import glob
-import io
 import os
+import warnings
 
-import joblib
 import numpy as np
-import six
 
 from polychrom.hdf5_format import load_URI
 
 from . import hdf5_format
 
 
 def load(filename):
-    """Universal load function for any type of data file It always returns just XYZ
-    positions - use fetch_block or hdf5_format.load_URI for loading the whole metadata
-
-    Accepted file types
-    -------------------
-
-    New-style URIs (HDF5 based storage)
-
-    Text files in openmm-polymer format
-    joblib files in openmm-polymer format
+    """
+    A function to load a single conformation from a URI. Deprecated.
+    Use load_URI from hdf5_format instead.
 
     Parameters
     ----------
@@ -59,37 +49,18 @@ def load(filename):
         filename to load or a URI
 
     """
+    warnings.warn("polymerutils.load is deprecated. Use hdf5_format.load_URI instead.", DeprecationWarning)
+
     if "::" in filename:
         return hdf5_format.load_URI(filename)["pos"]
 
-    if not os.path.exists(filename):
-        raise IOError("File not found :( \n %s" % filename)
-
-    try:  # loading from a joblib file here
-        return dict(joblib.load(filename)).pop("data")
-    except Exception:  # checking for a text file
-        data_file = open(filename)
-        line0 = data_file.readline()
-        try:
-            N = int(line0)
-        except (ValueError, UnicodeDecodeError):
-            raise TypeError("Could not read the file. Not text or joblib.")
-        data = [list(map(float, i.split())) for i in data_file.readlines()]
-
-        if len(data) != N:
-            raise ValueError("N does not correspond to the number of lines!")
-        return np.array(data)
+    raise ValueError("Only URIs are supported in this version of polychrom")
 
 
 def fetch_block(folder, ind, full_output=False):
     """
-    A more generic function to fetch block number "ind" from a trajectory in a folder
-
-
-    This function is useful both if you want to load both "old style" trajectories (block1.dat),
-    and "new style" trajectories ("blocks_1-50.h5")
-
-    It will be used in files "show"
+    A function to fetch a single block from a folder with a new-style trajectory.
+    Old-style trajectores are deprecated.
 
     Parameters
     ----------
@@ -109,12 +80,14 @@ def fetch_block(folder, ind, full_output=False):
 
         if full_output==True, then dict with data and metadata; XYZ is under key "pos"
     """
+    warnings.warn(
+        "fetch_block is deprecated. Use hdf5_format.list_uris followed by hdf5_format.load_URI instead.",
+        DeprecationWarning,
+    )
+
     blocksh5 = glob.glob(os.path.join(folder, "blocks*.h5"))
-    blocksdat = glob.glob(os.path.join(folder, "block*.dat"))
     ind = int(ind)
-    if (len(blocksh5) > 0) and (len(blocksdat) > 0):
-        raise ValueError("both .h5 and .dat files found in folder - exiting")
-    if (len(blocksh5) == 0) and (len(blocksdat) == 0):
+    if len(blocksh5) == 0:
         raise ValueError("no blocks found")
 
     if len(blocksh5) > 0:
@@ -129,45 +102,21 @@ def fetch_block(folder, ind, full_output=False):
         pos = exists.index(True)
         block = load_URI(blocksh5[pos] + f"::{ind}")
         if not full_output:
-            block = block["pos"]
+            return block["pos"]
+        return block
 
-    if len(blocksdat) > 0:
-        block = load(os.path.join(folder, f"block{ind}.dat"))
-    return block
+    raise ValueError(f"Cannot find the block {ind} in the folder {folder}")
 
 
 def save(data, filename, mode="txt", pdbGroups=None):
     """
-    Basically unchanged polymerutils.save function from openmm-polymer
-
-    It can save into txt or joblib formats used by old openmm-polymer
-
-    It is also very useful for saving files to PDB format to make them compatible
-    with nglview, pymol_show and others
+    A legacy function, currently only kept for compatibility with PDB saving that is rarely used.
     """
-    data = np.asarray(data, dtype=np.float32)
+    warnings.warn("polymerutils.save is deprecated. Will be moved to legacy", DeprecationWarning)
 
-    if mode.lower() == "joblib":
-        joblib.dump({"data": data}, filename=filename, compress=9)
-        return
-
-    if mode.lower() == "txt":
-        lines = [str(len(data)) + "\n"]
-
-        for particle in data:
-            lines.append("{0:.3f} {1:.3f} {2:.3f}\n".format(*particle))
-        if filename is None:
-            return lines
-
-        elif isinstance(filename, six.string_types):
-            with open(filename, "w") as myfile:
-                myfile.writelines(lines)
-        elif hasattr(filename, "writelines"):
-            filename.writelines(lines)
-        else:
-            raise ValueError("Not sure what to do with filename {0}".format(filename))
+    data = np.asarray(data, dtype=np.float32)
 
-    elif mode == "pdb":
+    if mode == "pdb":
         data = data - np.minimum(np.min(data, axis=0), np.zeros(3, float) - 100)[None, :]
         retret = ""
 
@@ -212,13 +161,10 @@ def add(st, n):
                 filename.write("C {0} {1} {2}".format(*i))
 
     else:
-        raise ValueError("Unknown mode : %s, use h5dict, joblib, txt or pdb" % mode)
+        raise ValueError(f"Unknown mode {mode}. Only 'pdb' and 'pyxyz' are supported.")
 
 
 def rotation_matrix(rotate):
-    """Calculates rotation matrix based on three rotation angles"""
-    tx, ty, tz = rotate
-    Rx = np.array([[1, 0, 0], [0, np.cos(tx), -np.sin(tx)], [0, np.sin(tx), np.cos(tx)]])
-    Ry = np.array([[np.cos(ty), 0, -np.sin(ty)], [0, 1, 0], [np.sin(ty), 0, np.cos(ty)]])
-    Rz = np.array([[np.cos(tz), -np.sin(tz), 0], [np.sin(tz), np.cos(tz), 0], [0, 0, 1]])
-    return np.dot(Rx, np.dot(Ry, Rz))
+    warnings.warn("rotation_matrix will be moved to polymer_analyses", DeprecationWarning)
+    from polychrom.polymer_analyses import rotation_matrix as rm
+    return rm(rotate)

pyproject.toml

@@ -1,9 +1,59 @@
+[build-system]
+requires = ["setuptools>=45", "wheel", "cython>=0.29", "numpy>=1.9"]
+build-backend = "setuptools.build_meta"
+
+[project]
+name = "polychrom"
+version = "0.1.1"
+description = "A library for polymer simulations and their analyses."
+readme = "README.md"
+requires-python = ">=3.7"
+license = "MIT"
+authors = [
+    {name = "Mirny Lab", email = "espresso@mit.edu"}
+]
+keywords = ["genomics", "polymer", "Hi-C", "molecular dynamics", "chromosomes"]
+classifiers = [
+    "Development Status :: 4 - Beta",
+    "Intended Audience :: Science/Research",
+    "Programming Language :: Python :: 3",
+    "Programming Language :: Python :: 3.7",
+    "Programming Language :: Python :: 3.8",
+    "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",
+]
+
+dependencies = [
+    "cython>=0.29",
+    "numpy>=1.9",
+    "scipy>=0.16",
+    "h5py>=2.5",
+    "pandas>=0.19",
+]
+
+[project.optional-dependencies]
+dev = [
+    "pytest",
+    "black",
+    "isort",
+]
+
+[project.urls]
+Homepage = "https://github.com/open2c/polychrom"
+Repository = "https://github.com/open2c/polychrom"
+
+[project.scripts]
+polychrom = "polychrom.cli:cli"
+
+[tool.setuptools.packages.find]
+include = ["polychrom*"]
+exclude = ["utilities*", "build*", "docs*", "examples*"]
+
 [tool.black]
 line-length = 120
 
 [tool.isort]
 profile = "black"
-
-[build-system]
-requires = ["setuptools", "setuptools-scm", "cython"]
-build-backend = "setuptools.build_meta"

requirements.txt

@@ -1,9 +1,6 @@
-six
 cython
 numpy>=1.9
 scipy>=0.16
 h5py>=2.5
 pandas>=0.19
-joblib
-pyknotid
 pytest