Repository: wangyiqiu/dbscan-python (language: python)

Issue:
pypi module not working
hello, the compiled module works like a charm but the module from pypi fails to load the DBSCAN method with error:
cannot import name 'DBSCAN' from 'dbscan' (unknown location)

awesome work btw
Cannot import in Conda environment
In my conda environment:
`>>> import dbscan
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/zpyang/anaconda3/envs/visnet/lib/python3.9/site-packages/dbscan/__init__.py", line 22, in <module>
    from dbscan.DBSCAN import DBSCAN
ModuleNotFoundError: No module named 'dbscan.DBSCAN'
`
but it is successfully imported while outside the conda environment. Not sure why that's the case.


Discussion / hints:
Thanks for your feedback. The module from pypi is kind of experimental at the moment, since we have only pre-compiled the library for Linux on x86_64. But we plan to extend it to support more computing environments in the near future.

May I know your computing environment? (e.g. architecture and operating system)
Thank you. 
it is arch linux on x86 using python 3.9.6

it looks like it fails to load the module but the path is in fact correct:
```
$ python -c "import dbscan"
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "/home/nik/.local/lib/python3.9/site-packages/dbscan/__init__.py", line 22, in <module>
    from dbscan.DBSCAN import DBSCAN
ModuleNotFoundError: No module named 'dbscan.DBSCAN'

$ ls /home/nik/.local/lib/python3.9/site-packages/dbscan
DBSCAN.cpython-38-x86_64-linux-gnu.so  example.py  __init__.py  __pycache__
```

if I change line 22 from `from dbscan.DBSCAN import DBSCAN` to `from dbscan import DBSCAN` the error changes to a circular import so my guess is that `__init__.py` tries to import the precompiled DBSCAN module and this tries to import `__init__.py`
Hi, thanks for the feedback. We have not tested the software in a Conda environment before. Could you let us know the system environment and Conda version that you are using?
I am using Ubuntu 20.04, conda 4.10.3, and python 3.9.7. Installed in dbscan using pip in conda environment. Hope this is useful. I also tried to build from source in conda environment but that didn't work as well. 

Fix patch:
diff --git a/README.md b/README.md
index acf4cda..9d9b99e 100644
--- a/README.md
+++ b/README.md
@@ -28,16 +28,16 @@ make -j # this will take a while
 
 The `<data-file>` can be any CSV-like point data file, where each line contains a data point -- see an example [here](https://github.com/wangyiqiu/hdbscan/blob/main/example-data.csv). The data file can be either with or without header. The cluster output `clusters.txt` will contain a cluster ID on each line (other than the first-line header), giving a cluster assignment in the same ordering as the input file. A noise point will have a cluster ID of `-1`.
 
-## Option 2: Use the Python binding (experimental)
+## Option 2: Use the Python binding
 
-We are developing a Python wrapper using Cython. Currently it is still experimental, as it is only tested to work on a fresh copy of Ubuntu 20.04 (x86_64 with Python 3.8+). There are two ways to install it:
+There are two ways to install it:
 
-* Compile it yourself: First install dependencies ``pip3 install -r src/requirements.txt`` and ``sudo apt install libpython3-dev``. Navigate to ``src/``, and run the Make script ``./make.sh``, The compilation will take a few minutes, and generate a ``.so`` library containing the ``DBSCAN`` module.
-* ***OR*** Install it using PyPI: ``pip3 install --user dbscan`` (the latest verion is 0.0.9)
+* Compile it yourself: First install dependencies ``pip3 install -r src/requirements.txt`` and ``sudo apt install libpython3-dev``. Run ``python3 setup.py build --inplace``, The compilation will take a few minutes, and generate a ``.so`` library containing the ``DBSCAN`` module.
+* ***OR*** Install it using PyPI: ``pip3 install --user dbscan`` (the latest version is 0.0.9)
 
 An example for using the Python module is provided in ``src/example.py``. If the dependencies above are installed, simply run ``python3 example.py`` from ``src/`` to reproduce the plots above.
 
-We plan to support more Python versions and environments, and we welcome any contributions that achieve these.
+To create a wheel that is supported universally across many Python versions for your given OS, run ``python setup.py bdist_wheel --py-limited-api=cp{YOUR_PYTHON_VERSION}`` in an environment containing the oldest numpy version available for the version of Python that you are compiling for. For example, for Python 3.8, use numpy 1.17 to compile the wheel. Then, the wheel will work on all Python and numpy versions that are newer that that for your given OS. This is done automatically when installing via pip. When creating the wheel, it is best for the user to correctly specify the version of Python they are using as the argument. It is done automatically, but the process messes with the internals of the wheel package which may not be stable.
 
 ### Python API
 
diff --git a/executable/main.cpp b/executable/main.cpp
index cdf8d6d..74c9f76 100644
--- a/executable/main.cpp
+++ b/executable/main.cpp
@@ -1,3 +1,5 @@
+#include <iostream>
+
 #include "dbscan/point.h"
 #include "dbscan/geometryIO.h"
 #include "dbscan/pbbs/parallel.h"
@@ -22,8 +24,8 @@ int main(int argc, char* argv[]) {
   int* cluster = new int[PIn.n / dim];
   double* data = PIn.A;
 
-  auto clusterer = Wrapper::Caller(data, dim, PIn.n / dim);
-  clusterer.computeDBSCAN(p_epsilon, p_minpts, coreFlag, cluster);
+  if (DBSCAN(data, dim, PIn.n / dim, p_epsilon, p_minpts, coreFlag, cluster))
+    cout << "Error: dimension >20 is not supported." << endl;
 
   if (oFile != NULL) {
     writeArrayToFile("cluster-id", cluster, PIn.n / dim, oFile);
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000..d5492a2
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,2 @@
+[build-system]
+requires = ["setuptools", "wheel", "oldest-supported-numpy"]
diff --git a/setup.py b/setup.py
index 2f66375..b3a838d 100644
--- a/setup.py
+++ b/setup.py
@@ -1,6 +1,23 @@
 import setuptools
 from setuptools.extension import Extension
+import numpy
+import sys
 
+# Override wheel's default ABI tag.
+try:
+    import wheel.bdist_wheel
+
+    current_tag = f"cp{sys.version_info.major}{sys.version_info.minor}"
+    old_initialize_options = wheel.bdist_wheel.bdist_wheel.initialize_options
+    def initialize_options(self):
+        old_initialize_options(self)
+        self.py_limited_api = current_tag
+    wheel.bdist_wheel.bdist_wheel.initialize_options = initialize_options
+except:
+    # Give up if it doesn't work. Not a big deal.
+    pass
+
+# Read README.md and set it as the description
 with open("README.md", "r", encoding="utf-8") as fh:
     long_description = fh.read()
 
@@ -15,19 +32,30 @@
     keywords='cluster clustering density dbscan',
     url="https://github.com/wangyiqiu/dbscan-python",
     license='MIT',
-    packages=[''],#setuptools.find_packages(),
-    package_dir={'': '.'},
-    package_data={'': ['dbscan/DBSCAN.cpython-38-x86_64-linux-gnu.so']},
+    packages=setuptools.find_packages(where='src'),
+    package_dir={'': 'src'},
+    ext_modules=[Extension(
+        "dbscan",
+        ["src/dbscanmodule.cpp", "src/Caller.cpp"],
+        language = 'c++',
+        extra_compile_args=["-Isrc", "-std=c++17", "-pthread", "-g", "-O3", "-fPIC"],
+        include_dirs=[numpy.get_include()],
+        py_limited_api=True,
+        define_macros=[('Py_LIMITED_API', '0x03020000'), ('NPY_NO_DEPRECATED_API', 'NPY_1_7_API_VERSION')]
+    )],
     classifiers=[
         'Development Status :: 2 - Pre-Alpha',
         'Intended Audience :: Science/Research',
         'Intended Audience :: Developers',
         "License :: OSI Approved :: MIT License",
         'Programming Language :: C++',
-        'Programming Language :: Python :: 3.8',
+        f'Programming Language :: Python :: 3',
         'Topic :: Software Development',
         'Topic :: Scientific/Engineering',
         "Operating System :: POSIX :: Linux",
     ],
-    python_requires='>=3.8',
+    python_requires=f'>={sys.version_info.major}.{sys.version_info.minor},<4',
+    install_requires=[
+       f'numpy>={numpy.__version__},<2'
+    ],
 )
diff --git a/src/Caller.cpp b/src/Caller.cpp
index fda06bc..f12a0dd 100644
--- a/src/Caller.cpp
+++ b/src/Caller.cpp
@@ -1,260 +1,34 @@
-#include <iostream>
 #include "Caller.h"
 
-#include "dbscan/point.h"
-#include "dbscan/shared.h"
-#include "dbscan/grid.h"
-#include "dbscan/coreBccp.h"
-#include "dbscan/pbbs/gettime.h"
+#include "dbscan/algo.h"
 #include "dbscan/pbbs/parallel.h"
-#include "dbscan/pbbs/sampleSort.h"
-#include "dbscan/pbbs/unionFind.h"
-
-// #define VERBOSE
-
-namespace Wrapper {
-
-  // Default constructor
-  Caller::Caller () {}
-
-  // Overloaded constructor
-  Caller::Caller (double* PF, intT dim, intT n) {
-    this->PF = PF;
-    this->n = n;
-    this->dim = dim;
-  }
-
-  // Destructor
-  Caller::~Caller () {}
-
-  template<int dim>
-  intT* computeDBSCANInternal(point<dim>* PRead, intT n, double epsilon, intT minPts, bool* coreFlagOut, intT* coreFlag, intT* cluster) {
-    typedef point<dim> pointT;
-    typedef grid<dim, pointT> gridT;
-    typedef cell<dim, pointT> cellT;
-
-#ifdef VERBOSE
-    cout << "Input: " << n << " points, dimension " << dim << endl;
-    printScheduler();
-    timing tt; tt.start();
-    timing t0; t0.start();
-#endif
-
-    floatT epsSqr = epsilon*epsilon;
-    pointT pMin = pMinParallel(PRead, n);
-
-    auto P = newA(pointT, n);
-    // parallel_for(0, n, [&](intT i){P[i] = PRead[i];});
-    auto G = new gridT(n+1, pMin, epsilon/sqrt(dim));
-    auto I = newA(intT, n);
-    G->insertParallel(PRead, P, n, I);
-#ifdef VERBOSE
-    cout << "num-cell = " << G->numCell() << endl;
-    cout << "compute-grid = " << t0.next() << endl;
-#endif
-    //mark core
-    parallel_for(0, n, [&](intT i) {coreFlag[i] = -1;});
-
-    auto isCore = [&](pointT *p) {
-		    coreFlag[p-P] = 1;
-		    return false;
-		  };
-
-    parallel_for(0, G->numCell(), [&](intT i) {
-	cellT* c = G->getCell(i);
-	if (c->size() >= minPts) c->pointMap(isCore);
-      });
-
-    parallel_for(0, n, [&](intT i) {
-			 if (coreFlag[i] < 0) {
-			   intT count = 0;
-			   auto isCore = [&] (pointT *p) {
-					   if(count >= minPts) return true;
-					   if(p->distSqr(P[i]) <= epsSqr) {//todo sqrt opt
-					     count ++;}
-					   return false;};
-			   G->nghPointMap(P[i].coordinate(), isCore);
-			   if (count >= minPts) coreFlag[i] = 1;
-			   else coreFlag[i] = 0;
-			 }
-		       });
-#ifdef VERBOSE
-    cout << "mark-core-time = " << t0.next() << endl;
-#endif
-    //cluster core
-    auto ccFlag = newA(intT, G->numCell());
-    parallel_for(0, G->numCell(), [&](intT i) {
-	auto ci = G->getCell(i);
-	ccFlag[i] = 0;
-	auto hasCore = [&](pointT *p) {
-			 if (coreFlag[p-P]) {
-			   ccFlag[i] = 1;
-			   return true;
-			 }
-			 return false;
-		       };
-	ci->pointMap(hasCore);
-      });
-
-    typedef kdTree<dim, pointT> treeT;
-    typedef kdNode<dim, pointT> nodeT;
-    typedef typename nodeT::bcp bcpT;
-    auto trees = newA(treeT*, G->numCell());
-    parallel_for(0, G->numCell(), [&](intT i) {trees[i] = NULL;});
-
-    // auto degCmp = [&](intT i, intT j) {
-    //                 return G->getCell(i)->size() < G->getCell(j)->size();
-    //               };
-    // auto ordering = newA(intT, G->numCell());
-    // par_for(intT i=0; i<G->numCell(); ++i) ordering[i] = i;
-    //sampleSort(ordering, G->numCell(), degCmp);
-
-    auto uf = unionFind(G->numCell());
-
-    floatT bcpTotalTime = 0; timing t1;
-    parallel_for(0, G->numCell(), [&](intT i) {
-	if (ccFlag[i]) {
-	  auto ti = trees[i];
-	  auto procTj = [&](cellT* cj) {
-			  intT j = cj - G->getCell(0);
-			  auto tj = trees[j];
-			  if (j < i && ccFlag[j] &&
-			      uf.find(i) != uf.find(j)) {
-			    if(hasEdge<cellT, treeT, pointT>(i, j, coreFlag, P, epsilon, G->getCell(0), trees)) {
-			      uf.link(i, j);
-			    }
-			  }
-			  return false;
-			};
-	  //G->nghCellMap(G->getCell(ordering[i]), procTj);
-	  G->nghCellMap(G->getCell(i), procTj);
-	}
-      });
-
-    parallel_for(0, G->numCell(), [&](intT i) {
-	if (trees[i]) delete trees[i];
-      });
-
-    parallel_for(0, n, [&](intT i) {cluster[i] = -1;});
-
-    parallel_for(0, G->numCell(), [&](intT i) {
-	auto cid = G->getCell(uf.find(i))->getItem() - P;//id of first point
-	auto clusterCore = [&](pointT* p){
-			     if (coreFlag[p - P])
-			       cluster[p - P] = cid;
-			     return false;
-			   };
-	G->getCell(i)->pointMap(clusterCore);
-      });
-#ifdef VERBOSE
-    cout << "cluster-core-time = " << t0.next() << endl;
-#endif
-    //cluster border to closest core point
-    parallel_for(0, n, [&](intT i) {
-			 if (!coreFlag[i]) {
-			   intT cid = -1;
-			   floatT cDistSqr = floatMax();
-			   auto closestCore = [&] (pointT* p) {
-						if (coreFlag[p-P]) {
-						  auto dist = p->distSqr(P[i]);
-						  if (dist <= epsSqr && dist < cDistSqr) {
-						    cDistSqr = dist;
-						    cid = cluster[p-P];}
-						}
-						return false;};
-			   G->nghPointMap(P[i].coordinate(), closestCore);
-			   cluster[i] = cid;
-			 }
-		       });
-#ifdef VERBOSE
-    cout << "cluster-border-time = " << t0.next() << endl;
-    cout << ">> total-clustering-time = " << tt.next() << endl;
-#endif
-    uf.del();
-    free(ccFlag);
-    free(trees);
-    delete G;
-
-    //improving cluster representation
-    auto cluster2 = newA(intT, n);
-    auto flag = newA(intT, n+1);
-    parallel_for(0, n, [&](intT i){cluster2[i] = cluster[i];});
-    sampleSort(cluster, n, std::less<intT>());
-
-    flag[0] = 1;
-    parallel_for(1, n, [&](intT i){
-			 if (cluster[i] != cluster[i-1])
-			   flag[i] = 1;
-			 else
-			   flag[i] = 0;
-		       });
-    flag[n] = sequence::prefixSum(flag, 0, n);
-
-    typedef Table<hashSimplePair,intT> tableT;
-    auto T = new tableT(n, hashSimplePair());
-    parallel_for(0, n, [&] (intT i) {
-			 if (flag[i] != flag[i+1]) {
-			   T->insert(make_pair(cluster[i], flag[i]));
-			 }
-		       });
-
-    if(T->find(-1).second < 0) {
-      parallel_for(0, n, [&](intT i){
-			   cluster2[i] = T->find(cluster2[i]).second;
-			 });
-    } else {
-      parallel_for(0, n, [&](intT i){
-			   if (cluster2[i] > 0)
-			     cluster2[i] = T->find(cluster2[i]).second-1;
-			 });
-    }
-
-    //restoring order
-    parallel_for(0, n, [&](intT i){
-			 cluster[I[i]] = cluster2[i];
-		       });
-    parallel_for(0, n, [&](intT i){
-			 coreFlagOut[I[i]] = coreFlag[i];
-		       });
-
-    free(I);
-    free(cluster2);
-    free(flag);
-    T->del(); // Required to clean-up T's internals
-    delete T;
-    free(P);
-#ifdef VERBOSE
-    cout << "output-time = " << tt.stop() << endl;
-#endif
-    return cluster;
+#include "dbscan/pbbs/utils.h"
+
+int DBSCAN(floatT* PF, intT dim, intT n, double epsilon, intT minPts, bool* coreFlag, intT* labels) {
+  auto coreFlag2 = newA(intT, n);
+  if (dim == 2) {DBSCAN<2>(PF, n, epsilon, minPts, coreFlag, coreFlag2, labels);}
+  else if (dim == 3) {DBSCAN<3>(PF, n, epsilon, minPts, coreFlag, coreFlag2, labels);}
+  else if (dim == 4) {DBSCAN<4>(PF, n, epsilon, minPts, coreFlag, coreFlag2, labels);}
+  else if (dim == 5) {DBSCAN<5>(PF, n, epsilon, minPts, coreFlag, coreFlag2, labels);}
+  else if (dim == 6) {DBSCAN<6>(PF, n, epsilon, minPts, coreFlag, coreFlag2, labels);}
+  else if (dim == 7) {DBSCAN<7>(PF, n, epsilon, minPts, coreFlag, coreFlag2, labels);}
+  else if (dim == 8) {DBSCAN<8>(PF, n, epsilon, minPts, coreFlag, coreFlag2, labels);}
+  else if (dim == 9) {DBSCAN<9>(PF, n, epsilon, minPts, coreFlag, coreFlag2, labels);}
+  else if (dim == 10) {DBSCAN<10>(PF, n, epsilon, minPts, coreFlag, coreFlag2, labels);}
+  else if (dim == 11) {DBSCAN<11>(PF, n, epsilon, minPts, coreFlag, coreFlag2, labels);}
+  else if (dim == 12) {DBSCAN<12>(PF, n, epsilon, minPts, coreFlag, coreFlag2, labels);}
+  else if (dim == 13) {DBSCAN<13>(PF, n, epsilon, minPts, coreFlag, coreFlag2, labels);}
+  else if (dim == 14) {DBSCAN<14>(PF, n, epsilon, minPts, coreFlag, coreFlag2, labels);}
+  else if (dim == 15) {DBSCAN<15>(PF, n, epsilon, minPts, coreFlag, coreFlag2, labels);}
+  else if (dim == 16) {DBSCAN<16>(PF, n, epsilon, minPts, coreFlag, coreFlag2, labels);}
+  else if (dim == 17) {DBSCAN<17>(PF, n, epsilon, minPts, coreFlag, coreFlag2, labels);}
+  else if (dim == 18) {DBSCAN<18>(PF, n, epsilon, minPts, coreFlag, coreFlag2, labels);}
+  else if (dim == 19) {DBSCAN<19>(PF, n, epsilon, minPts, coreFlag, coreFlag2, labels);}
+  else if (dim == 20) {DBSCAN<20>(PF, n, epsilon, minPts, coreFlag, coreFlag2, labels);}
+  else {
+    return 1;
   }
-
-  intT* Caller::computeDBSCAN(double epsilon, intT minPts, bool* coreFlag, intT* labels) {
-    auto coreFlag2 = newA(intT, n);
-    if (this->dim == 2) {computeDBSCANInternal<2>((point<2>*) PF, this->n, epsilon, minPts, coreFlag, coreFlag2, labels);}
-    else if (this->dim == 3) {computeDBSCANInternal<3>((point<3>*) PF, this->n, epsilon, minPts, coreFlag, coreFlag2, labels);}
-    else if (this->dim == 4) {computeDBSCANInternal<4>((point<4>*) PF, this->n, epsilon, minPts, coreFlag, coreFlag2, labels);}
-    else if (this->dim == 5) {computeDBSCANInternal<5>((point<5>*) PF, this->n, epsilon, minPts, coreFlag, coreFlag2, labels);}
-    else if (this->dim == 6) {computeDBSCANInternal<6>((point<6>*) PF, this->n, epsilon, minPts, coreFlag, coreFlag2, labels);}
-    else if (this->dim == 7) {computeDBSCANInternal<7>((point<7>*) PF, this->n, epsilon, minPts, coreFlag, coreFlag2, labels);}
-    else if (this->dim == 8) {computeDBSCANInternal<8>((point<8>*) PF, this->n, epsilon, minPts, coreFlag, coreFlag2, labels);}
-    else if (this->dim == 9) {computeDBSCANInternal<9>((point<9>*) PF, this->n, epsilon, minPts, coreFlag, coreFlag2, labels);}
-    else if (this->dim == 10) {computeDBSCANInternal<10>((point<10>*) PF, this->n, epsilon, minPts, coreFlag, coreFlag2, labels);}
-    else if (this->dim == 11) {computeDBSCANInternal<11>((point<11>*) PF, this->n, epsilon, minPts, coreFlag, coreFlag2, labels);}
-    else if (this->dim == 12) {computeDBSCANInternal<12>((point<12>*) PF, this->n, epsilon, minPts, coreFlag, coreFlag2, labels);}
-    else if (this->dim == 13) {computeDBSCANInternal<13>((point<13>*) PF, this->n, epsilon, minPts, coreFlag, coreFlag2, labels);}
-    else if (this->dim == 14) {computeDBSCANInternal<14>((point<14>*) PF, this->n, epsilon, minPts, coreFlag, coreFlag2, labels);}
-    else if (this->dim == 15) {computeDBSCANInternal<15>((point<15>*) PF, this->n, epsilon, minPts, coreFlag, coreFlag2, labels);}
-    else if (this->dim == 16) {computeDBSCANInternal<16>((point<16>*) PF, this->n, epsilon, minPts, coreFlag, coreFlag2, labels);}
-    else if (this->dim == 17) {computeDBSCANInternal<17>((point<17>*) PF, this->n, epsilon, minPts, coreFlag, coreFlag2, labels);}
-    else if (this->dim == 18) {computeDBSCANInternal<18>((point<18>*) PF, this->n, epsilon, minPts, coreFlag, coreFlag2, labels);}
-    else if (this->dim == 19) {computeDBSCANInternal<19>((point<19>*) PF, this->n, epsilon, minPts, coreFlag, coreFlag2, labels);}
-    else if (this->dim == 20) {computeDBSCANInternal<20>((point<20>*) PF, this->n, epsilon, minPts, coreFlag, coreFlag2, labels);}
-    else {
-      cout << "Error: dimension >20 is not supported." << endl;
-    }
-    free(coreFlag2);
-    return labels;
-  }
-
+  free(coreFlag2);
+  return 0;
 }
+
diff --git a/src/Caller.h b/src/Caller.h
index 4a03195..1268077 100644
--- a/src/Caller.h
+++ b/src/Caller.h
@@ -1,16 +1,7 @@
 #pragma once
 
-#include "dbscan/pbbs/parallel.h"
+// intT* DBSCAN(floatT* PF, int dim, intT n, floatT epsilon, intT minPts, bool* coreFlag, intT* cluster);
 
-namespace Wrapper {
-  class Caller {
-  public:
-    intT n;
-    int dim;
-    floatT* PF;
-    Caller();
-    Caller(floatT* PF, int dim, intT n);
-    ~Caller();
-    intT* computeDBSCAN(floatT epsilon, intT minPts, bool* coreFlag, intT* cluster);
-  };
-}
+// replace names from "dbscan/pbbs/parallel.h" with actual types
+
+int DBSCAN(double* PF, int dim, int n, double epsilon, int minPts, bool* coreFlag, int* cluster);
diff --git a/src/Caller.pxd b/src/Caller.pxd
deleted file mode 100644
index 4523110..0000000
--- a/src/Caller.pxd
+++ /dev/null
@@ -1,11 +0,0 @@
-from libcpp cimport bool
-
-cdef extern from "Caller.cpp":
-    pass
-
-# Declare the class with cdef
-cdef extern from "Caller.h" namespace "Wrapper":
-    cdef cppclass Caller:
-        Caller() except +
-        Caller(double *, int, int) except +
-        int* computeDBSCAN(double, int, bool *, int *)
diff --git a/src/DBSCAN.pyx b/src/DBSCAN.pyx
deleted file mode 100644
index a6bf6d9..0000000
--- a/src/DBSCAN.pyx
+++ /dev/null
@@ -1,37 +0,0 @@
-# distutils: language = c++
-
-from libc.stdlib cimport malloc, free
-from libcpp cimport bool
-from Caller cimport Caller
-import numpy as np
-
-def DBSCAN(X, double eps=0.5, int min_samples=5):
-    cdef int n = X.shape[0]
-    dim = X.shape[1]
-
-    if dim <= 1:
-        print("Error: invalid input data dimensionality (has to >1).")
-        return np.array([],dtype=np.int32), np.array([],dtype=np.bool)
-
-    if n <= 0:
-        return np.array([],dtype=np.int32), np.array([],dtype=np.bool)
-
-    if n > 100000000:
-        print("Warning: large n, the program behavior might be undefined due to overflow.")
-
-    cdef double [:,:] X_view = X
-
-    my_instance = new Caller(&X_view[0,0], dim, n)
-
-    core_samples = np.empty(n, dtype=np.bool)
-    cdef bool [:] core_samples_view = core_samples
-
-    labels = np.empty(n, dtype=np.int32)
-    cdef int [:] labels_view = labels
-
-    try:
-        my_instance.computeDBSCAN(eps, min_samples, &core_samples_view[0], &labels_view[0])
-    finally:
-        del my_instance
-
-    return labels, core_samples
diff --git a/src/__init__.py b/src/__init__.py
deleted file mode 100644
index 3889d58..0000000
--- a/src/__init__.py
+++ /dev/null
@@ -1,22 +0,0 @@
-platformOk = True
-
-import platform
-
-if platform.system() != 'Linux':
-    print("Detected an operating system other than Linux, the DBSCAN package only works with 64-bit Linux at the moment.")
-    platformOk = False
-
-import sys
-
-if sys.version_info[0] < 3:
-    print("Detected Python " + str(sys.version_info[0]) + ", the DBSCAN package requires Python 3.8+.")
-    platformOk = False
-
-if sys.version_info[0] == 3 and sys.version_info[1] < 8:
-    print("Detected Python " + str(sys.version_info[0]) + "." + str(sys.version_info[1]) + ", the DBSCAN package requires Python 3.8+")
-    platformOk = False
-
-if not platformOk:
-    print("Please refer to https://github.com/wangyiqiu/dbscan-python for more specific system requirements. Please also make sure the latest version of the package (0.0.9) is installed.")
-else:
-    from dbscan.DBSCAN import DBSCAN
diff --git a/src/clean.sh b/src/clean.sh
deleted file mode 100755
index a9525ce..0000000
--- a/src/clean.sh
+++ /dev/null
@@ -1,7 +0,0 @@
-#!/bin/sh
-rm *.o
-rm *.so
-rm *~
-rm DBSCAN.cpp
-rm -r build
-rm *.html
diff --git a/src/dbscan/algo.h b/src/dbscan/algo.h
new file mode 100644
index 0000000..6699aed
--- /dev/null
+++ b/src/dbscan/algo.h
@@ -0,0 +1,218 @@
+#pragma once
+
+#include <iostream>
+#include "dbscan/point.h"
+#include "dbscan/shared.h"
+#include "dbscan/grid.h"
+#include "dbscan/coreBccp.h"
+#include "dbscan/pbbs/gettime.h"
+#include "dbscan/pbbs/parallel.h"
+#include "dbscan/pbbs/sampleSort.h"
+#include "dbscan/pbbs/unionFind.h"
+
+// #define VERBOSE
+
+template<int dim>
+intT* DBSCAN(floatT* PF, intT n, double epsilon, intT minPts, bool* coreFlagOut, intT* coreFlag, intT* cluster) {
+  typedef point<dim> pointT;
+  typedef grid<dim, pointT> gridT;
+  typedef cell<dim, pointT> cellT;
+
+  point<dim>* PRead = (point<dim>*)PF;
+
+#ifdef VERBOSE
+  cout << "Input: " << n << " points, dimension " << dim << endl;
+  printScheduler();
+  timing tt; tt.start();
+  timing t0; t0.start();
+#endif
+
+  floatT epsSqr = epsilon*epsilon;
+  pointT pMin = pMinParallel(PRead, n);
+
+  auto P = newA(pointT, n);
+  // parallel_for(0, n, [&](intT i){P[i] = PRead[i];});
+  auto G = new gridT(n+1, pMin, epsilon/sqrt(dim));
+  auto I = newA(intT, n);
+  G->insertParallel(PRead, P, n, I);
+#ifdef VERBOSE
+  cout << "num-cell = " << G->numCell() << endl;
+  cout << "compute-grid = " << t0.next() << endl;
+#endif
+  //mark core
+  parallel_for(0, n, [&](intT i) {coreFlag[i] = -1;});
+
+  auto isCore = [&](pointT *p) {
+                  coreFlag[p-P] = 1;
+                  return false;
+                };
+
+  parallel_for(0, G->numCell(), [&](intT i) {
+      cellT* c = G->getCell(i);
+      if (c->size() >= minPts) c->pointMap(isCore);
+    });
+
+  parallel_for(0, n, [&](intT i) {
+                       if (coreFlag[i] < 0) {
+                         intT count = 0;
+                         auto isCore = [&] (pointT *p) {
+                                         if(count >= minPts) return true;
+                                         if(p->distSqr(P[i]) <= epsSqr) {//todo sqrt opt
+                                           count ++;}
+                                         return false;};
+                         G->nghPointMap(P[i].coordinate(), isCore);
+                         if (count >= minPts) coreFlag[i] = 1;
+                         else coreFlag[i] = 0;
+                       }
+                     });
+#ifdef VERBOSE
+  cout << "mark-core-time = " << t0.next() << endl;
+#endif
+  //cluster core
+  auto ccFlag = newA(intT, G->numCell());
+  parallel_for(0, G->numCell(), [&](intT i) {
+      auto ci = G->getCell(i);
+      ccFlag[i] = 0;
+      auto hasCore = [&](pointT *p) {
+                       if (coreFlag[p-P]) {
+                         ccFlag[i] = 1;
+                         return true;
+                       }
+                       return false;
+                     };
+      ci->pointMap(hasCore);
+    });
+
+  typedef kdTree<dim, pointT> treeT;
+  typedef kdNode<dim, pointT> nodeT;
+  typedef typename nodeT::bcp bcpT;
+  auto trees = newA(treeT*, G->numCell());
+  parallel_for(0, G->numCell(), [&](intT i) {trees[i] = NULL;});
+
+  // auto degCmp = [&](intT i, intT j) {
+  //                 return G->getCell(i)->size() < G->getCell(j)->size();
+  //               };
+  // auto ordering = newA(intT, G->numCell());
+  // par_for(intT i=0; i<G->numCell(); ++i) ordering[i] = i;
+  //sampleSort(ordering, G->numCell(), degCmp);
+
+  auto uf = unionFind(G->numCell());
+
+  floatT bcpTotalTime = 0; timing t1;
+  parallel_for(0, G->numCell(), [&](intT i) {
+      if (ccFlag[i]) {
+        auto ti = trees[i];
+        auto procTj = [&](cellT* cj) {
+                        intT j = cj - G->getCell(0);
+                        auto tj = trees[j];
+                        if (j < i && ccFlag[j] &&
+                            uf.find(i) != uf.find(j)) {
+                          if(hasEdge<cellT, treeT, pointT>(i, j, coreFlag, P, epsilon, G->getCell(0), trees)) {
+                            uf.link(i, j);
+                          }
+                        }
+                        return false;
+                      };
+        //G->nghCellMap(G->getCell(ordering[i]), procTj);
+        G->nghCellMap(G->getCell(i), procTj);
+      }
+    });
+
+  parallel_for(0, G->numCell(), [&](intT i) {
+      if (trees[i]) delete trees[i];
+    });
+
+  parallel_for(0, n, [&](intT i) {cluster[i] = -1;});
+
+  parallel_for(0, G->numCell(), [&](intT i) {
+      auto cid = G->getCell(uf.find(i))->getItem() - P;//id of first point
+      auto clusterCore = [&](pointT* p){
+                           if (coreFlag[p - P])
+                             cluster[p - P] = cid;
+                           return false;
+                         };
+      G->getCell(i)->pointMap(clusterCore);
+    });
+#ifdef VERBOSE
+  cout << "cluster-core-time = " << t0.next() << endl;
+#endif
+  //cluster border to closest core point
+  parallel_for(0, n, [&](intT i) {
+                       if (!coreFlag[i]) {
+                         intT cid = -1;
+                         floatT cDistSqr = floatMax();
+                         auto closestCore = [&] (pointT* p) {
+                                              if (coreFlag[p-P]) {
+                                                auto dist = p->distSqr(P[i]);
+                                                if (dist <= epsSqr && dist < cDistSqr) {
+                                                  cDistSqr = dist;
+                                                  cid = cluster[p-P];}
+                                              }
+                                              return false;};
+                         G->nghPointMap(P[i].coordinate(), closestCore);
+                         cluster[i] = cid;
+                       }
+                     });
+#ifdef VERBOSE
+  cout << "cluster-border-time = " << t0.next() << endl;
+  cout << ">> total-clustering-time = " << tt.next() << endl;
+#endif
+  uf.del();
+  free(ccFlag);
+  free(trees);
+  delete G;
+
+  //improving cluster representation
+  auto cluster2 = newA(intT, n);
+  auto flag = newA(intT, n+1);
+  parallel_for(0, n, [&](intT i){cluster2[i] = cluster[i];});
+  sampleSort(cluster, n, std::less<intT>());
+
+  flag[0] = 1;
+  parallel_for(1, n, [&](intT i){
+                       if (cluster[i] != cluster[i-1])
+                         flag[i] = 1;
+                       else
+                         flag[i] = 0;
+                     });
+  flag[n] = sequence::prefixSum(flag, 0, n);
+
+  typedef Table<hashSimplePair,intT> tableT;
+  auto T = new tableT(n, hashSimplePair());
+  parallel_for(0, n, [&] (intT i) {
+                       if (flag[i] != flag[i+1]) {
+                         T->insert(make_pair(cluster[i], flag[i]));
+                       }
+                     });
+
+  if(T->find(-1).second < 0) {
+    parallel_for(0, n, [&](intT i){
+                         cluster2[i] = T->find(cluster2[i]).second;
+                       });
+  } else {
+    parallel_for(0, n, [&](intT i){
+                         if (cluster2[i] > 0)
+                           cluster2[i] = T->find(cluster2[i]).second-1;
+                       });
+  }
+
+  //restoring order
+  parallel_for(0, n, [&](intT i){
+                       cluster[I[i]] = cluster2[i];
+                     });
+  parallel_for(0, n, [&](intT i){
+                       coreFlagOut[I[i]] = coreFlag[i];
+                     });
+
+  free(I);
+  free(cluster2);
+  free(flag);
+  T->del(); // Required to clean-up T's internals
+  delete T;
+  free(P);
+#ifdef VERBOSE
+  cout << "output-time = " << tt.stop() << endl;
+#endif
+  return cluster;
+}
+
diff --git a/src/dbscanmodule.cpp b/src/dbscanmodule.cpp
new file mode 100644
index 0000000..4a2940b
--- /dev/null
+++ b/src/dbscanmodule.cpp
@@ -0,0 +1,124 @@
+#include "Python.h"
+#include "numpy/arrayobject.h"
+#include "Caller.h"
+
+static PyObject* DBSCAN(PyObject* self, PyObject* args, PyObject *kwargs)
+{
+    PyObject *Xobj;
+    PyArrayObject *X = NULL;
+    double eps = 0.5;
+    int min_samples = 5;
+
+    static char *kwlist[] = {"X", "eps", "min_samples", NULL};
+
+    if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O|di", kwlist,
+                                     &Xobj, &eps, &min_samples))
+    {
+        return NULL;
+    }
+
+    // Check the number of dimensions and that we actually recieved an np.ndarray
+    X = (PyArrayObject*)PyArray_FROMANY(
+        Xobj,
+        NPY_DOUBLE,
+        2,
+        2,
+        NPY_ARRAY_CARRAY_RO
+    );
+
+    if (X == NULL)
+    {
+        return NULL;
+    }
+
+    npy_intp *dims = PyArray_SHAPE(X);
+
+    npy_intp n = dims[0];
+    npy_intp dim = dims[1];
+
+    if (dim <= 1)
+    {
+        PyErr_SetString(PyExc_ValueError, "DBSCAN: invalid input data dimensionality (has to >1)");
+        return NULL;
+    }
+
+    if (dim > 20)
+    {
+        PyErr_SetString(PyExc_ValueError, "DBSCAN: dimension >20 is not supported");
+        return NULL;
+    }
+
+    if (n > 100000000)
+    {
+        PyErr_WarnEx(PyExc_RuntimeWarning, "DBSCAN: large n, the program behavior might be undefined due to overflow", 1);
+    }
+
+    PyArrayObject* core_samples = (PyArrayObject*)PyArray_SimpleNew(1, &n, NPY_BOOL);
+    PyArrayObject* labels = (PyArrayObject*)PyArray_SimpleNew(1, &n, NPY_INT);
+
+    DBSCAN(
+        (double*)PyArray_DATA(X),
+        dim,
+        n,
+        eps,
+        min_samples,
+        (bool*)PyArray_DATA(core_samples),
+        (int*)PyArray_DATA(labels)
+    );
+
+    PyObject* ret = PyTuple_Pack(2, labels, core_samples);
+    Py_IncRef(ret);
+    return ret;
+}
+
+PyDoc_STRVAR(doc_DBSCAN,
+"DBSCAN(X, eps=0.5, min_samples=5)\n--\n\n\
+Run DBSCAN on a set of n samples of dimension dim with a minimum seperation\n\
+between the clusters (which must include at least min_samples) of eps. Points\n\
+that do not fit in any cluster are labeled as noise (-1).\n\
+\n\
+This function returns a tuple consisting of an int array of length n containing\n\
+the labels after clustering and a bool array of length n which differentiates\n\
+whether or not the sample is the core sample of its cluster.\n\
+\n\
+Parameters\n\
+----------\n\
+X : np.ndarray[tuple[n, dim], np.float64]\n\
+    2-D array representing the samples.\n\
+eps : float\n\
+    minimum seperation between the clusters.\n\
+min_samples : int\n\
+    minimum number of samples in the clusters.\n\
+\n\
+Returns\n\
+-------\n\
+labels : np.ndarray[tuple[n], np.int_]\n\
+    the labels after clustering\n\
+core_samples : np.ndarray[tuple[n], np.bool_]\n\
+    is each sample the core sample of its cluster\n\
+\n");
+
+static struct PyMethodDef methods[] = {
+    {"DBSCAN", (PyCFunction)(void*)(PyCFunctionWithKeywords) DBSCAN, METH_VARARGS | METH_KEYWORDS, doc_DBSCAN},
+    {NULL, NULL, 0, NULL}
+};
+
+typedef struct {
+    PyObject *DBSCAN;
+} dbscanModuleState;
+
+static struct PyModuleDef dbscanModule =
+{
+    PyModuleDef_HEAD_INIT,
+    "dbscan",
+    "",
+    sizeof(dbscanModuleState),
+    methods
+};
+
+PyMODINIT_FUNC
+PyInit_dbscan (void)
+{
+    import_array();
+    return PyModule_Create(&dbscanModule);
+}
diff --git a/src/example.py b/src/example.py
index fecc691..7480372 100644
--- a/src/example.py
+++ b/src/example.py
@@ -20,7 +20,7 @@
     core_samples_mask[db.core_sample_indices_] = True
     labels = db.labels_
 else:
-    from DBSCAN import DBSCAN
+    from dbscan import DBSCAN
     labels, core_samples_mask = DBSCAN(X, eps=0.3, min_samples=10)
 
 # #############################################################################
diff --git a/src/make.sh b/src/make.sh
deleted file mode 100755
index 2dc693e..0000000
--- a/src/make.sh
+++ /dev/null
@@ -1,5 +0,0 @@
-#!/bin/sh
-sh clean.sh
-cython -a --cplus DBSCAN.pyx
-gcc -std=c++17 -pthread -g -O3 -fPIC -I. -I/usr/include/python3.8 -c DBSCAN.cpp -o DBSCAN.cpython-38-x86_64-linux-gnu.o
-g++ -std=c++17 -pthread -shared -O3 DBSCAN.cpython-38-x86_64-linux-gnu.o -o DBSCAN.cpython-38-x86_64-linux-gnu.so


Test patch:
diff --git a/test/dbscan_test.cpp b/test/dbscan_test.cpp
index 875acdf..6265e5c 100644
--- a/test/dbscan_test.cpp
+++ b/test/dbscan_test.cpp
@@ -17,8 +17,7 @@ TEST(testDBSCAN, smallCluster1) {
   data[6] = 2.5; data[7] = 1.5;
   data[8] = 4; data[9] = 0;
 
-  auto clusterer = Wrapper::Caller(data.get(), dim, n);
-  clusterer.computeDBSCAN(1.42, 3, coreFlag.get(), cluster.get());
+  DBSCAN(data.get(), dim, n, 1.42, 3, coreFlag.get(), cluster.get());
 
   EXPECT_EQ(coreFlag[0], 0);
   EXPECT_EQ(coreFlag[1], 1);
@@ -46,8 +45,7 @@ TEST(testDBSCAN, smallCluster1_largeEps) {
   data[6] = 2.5; data[7] = 1.5;
   data[8] = 4; data[9] = 0;
 
-  auto clusterer = Wrapper::Caller(data.get(), dim, n);
-  clusterer.computeDBSCAN(5.7, 3, coreFlag.get(), cluster.get());
+  DBSCAN(data.get(), dim, n, 5.7, 3, coreFlag.get(), cluster.get());
 
   EXPECT_EQ(coreFlag[0], 1);
   EXPECT_EQ(coreFlag[1], 1);
@@ -75,8 +73,7 @@ TEST(testDBSCAN, smallCluster1_smallEps) {
   data[6] = 2.5; data[7] = 1.5;
   data[8] = 4; data[9] = 0;
 
-  auto clusterer = Wrapper::Caller(data.get(), dim, n);
-  clusterer.computeDBSCAN(0.7, 3, coreFlag.get(), cluster.get());
+  DBSCAN(data.get(), dim, n, 0.7, 3, coreFlag.get(), cluster.get());
 
   EXPECT_EQ(coreFlag[0], 0);
   EXPECT_EQ(coreFlag[1], 0);