custom clustering now takes position matrix instead of distance^2 matrix

pypolychord/_pypolychord.cpp

@@ -117,24 +117,24 @@ void dumper(int ndead, int nlive, int npars, double* live, double* dead, double*
 /* Callback to the cluster */
 static PyObject *python_cluster = NULL;
 
-void cluster(double* distance2_matrix, int* cluster_list, int n)
+void cluster(double* position_matrix, int* cluster_list, int m, int n)
 {
-    /* create a python version of distance2_matrix */
-    npy_intp shape[] = {n,n};            
-    PyObject *array_distance2_matrix = PyArray_SimpleNewFromData(2, shape, NPY_DOUBLE, distance2_matrix);
-    if (array_distance2_matrix ==NULL) throw PythonException();
-    PyArray_CLEARFLAGS(reinterpret_cast<PyArrayObject*>(array_distance2_matrix), NPY_ARRAY_WRITEABLE);
+    /* create a python version of position_matrix */
+    npy_intp shape[] = {m,n};            
+    PyObject *array_position_matrix = PyArray_SimpleNewFromData(2, shape, NPY_DOUBLE, position_matrix);
+    if (array_position_matrix ==NULL) throw PythonException();
+    PyArray_CLEARFLAGS(reinterpret_cast<PyArrayObject*>(array_position_matrix), NPY_ARRAY_WRITEABLE);
 
     /* create a python version of cluster_list */
     npy_intp shape1[] = {n};            
     PyObject *array_cluster_list = PyArray_SimpleNewFromData(1, shape1, NPY_INT, cluster_list);
-    if (array_cluster_list==NULL) {Py_DECREF(array_distance2_matrix); throw PythonException();}
+    if (array_cluster_list==NULL) {Py_DECREF(array_position_matrix); throw PythonException();}
 
     /* Compute cluster_list from the cluster */
-    PyObject_CallFunctionObjArgs(python_cluster,array_distance2_matrix,array_cluster_list,NULL);
+    PyObject_CallFunctionObjArgs(python_cluster,array_position_matrix,array_cluster_list,NULL);
 
     /* Garbage collect */
-    Py_DECREF(array_cluster_list); Py_DECREF(array_distance2_matrix);
+    Py_DECREF(array_cluster_list); Py_DECREF(array_position_matrix);
 }
 
 

pypolychord/polychord.py

@@ -12,8 +12,8 @@ def default_dumper(live, dead, logweights, logZ, logZerr):
     pass
 
 
-def default_cluster(distance2_matrix):
-    return np.zeros(distance2_matrix.shape[0],dtype=int)
+def default_cluster(position_matrix):
+    return np.zeros(position_matrix.shape[1],dtype=int)
 
 
 def run_polychord(loglikelihood, nDims, nDerived, settings,
@@ -112,8 +112,8 @@ def run_polychord(loglikelihood, nDims, nDerived, settings,
 
         Parameters
         ----------
-        distance2_matrix: numpy.array
-            squared distances between points. Shape (nPoints, nPoints)
+        position_matrix: numpy.array
+            positions of points. Shape (nDims, nPoints)
 
         Returns
         -------
@@ -197,8 +197,8 @@ def wrap_loglikelihood(theta, phi):
     def wrap_prior(cube, theta):
         theta[:] = prior(cube)
 
-    def wrap_cluster(distance2_matrix, cluster_list):
-        cluster_list[:] = cluster(distance2_matrix)
+    def wrap_cluster(position_matrix, cluster_list):
+        cluster_list[:] = cluster(position_matrix)
 
     # Run polychord from module library
     _pypolychord.run(wrap_loglikelihood,

run_pypolychord.py

@@ -41,8 +41,8 @@ def dumper(live, dead, logweights, logZ, logZerr):
 
 #| Optional cluster function allow user-defined clustering
 
-def cluster(distance2_matrix):
-    npoints = distance2_matrix.shape[0]
+def cluster(position_matrix):
+    npoints = position_matrix.shape[1]
     clusters = np.ones(npoints, dtype=int)
 
     # <do some clustering algorithm to assign clusters>

src/polychord/calculate.f90

@@ -91,22 +91,22 @@ end function calculate_posterior_point
     !! The final term can be written as a data_array^T data_array, and the first
     !! two are easy to write. We can therefore calculate this in two lines with
     !! instrisic functions
-    function calculate_similarity_matrix(data_array) result(similarity_matrix)
+    function calculate_distance2_matrix(data_array) result(distance2_matrix)
 
         real(dp), intent(in), dimension(:,:) :: data_array
 
-        real(dp), dimension(size(data_array,2),size(data_array,2)) :: similarity_matrix
+        real(dp), dimension(size(data_array,2),size(data_array,2)) ::distance2_matrix 
 
         integer :: i
 
 
-        similarity_matrix = spread( &
+        distance2_matrix = spread( &
             [ ( dot_product(data_array(:,i),data_array(:,i)), i=1,size(data_array,2) ) ], &
             dim=2,ncopies=size(data_array,2) )
 
-        similarity_matrix = similarity_matrix + transpose(similarity_matrix) - 2d0 * matmul( transpose(data_array),data_array )
+        distance2_matrix = distance2_matrix + transpose(distance2_matrix) - 2d0 * matmul( transpose(data_array),data_array )
 
-    end function calculate_similarity_matrix
+    end function calculate_distance2_matrix
 
 
 

src/polychord/clustering.f90

@@ -254,7 +254,7 @@ module cluster_module
     function do_clustering(settings,RTI,cluster,sub_dimensions)
         use settings_module,   only: program_settings
         use run_time_module,   only: run_time_info,add_cluster
-        use calculate_module,  only: calculate_similarity_matrix
+        use calculate_module,  only: calculate_distance2_matrix
         use KNN_clustering,    only: NN_clustering
         implicit none
 
@@ -266,10 +266,10 @@ function do_clustering(settings,RTI,cluster,sub_dimensions)
         integer,dimension(:),optional,intent(in) :: sub_dimensions
 
         interface
-            function cluster(distance2_matrix) result(cluster_list)
+            function cluster(position_matrix) result(cluster_list)
                 import :: dp
-                real(dp), intent(in), dimension(:,:) :: distance2_matrix
-                integer, dimension(size(distance2_matrix,1)) :: cluster_list
+                real(dp), intent(in), dimension(:,:) :: position_matrix
+                integer, dimension(size(position_matrix,1)) :: cluster_list
             end function
         end interface
 
@@ -278,14 +278,16 @@ function cluster(distance2_matrix) result(cluster_list)
 
 
         ! Similarity matrix
-        real(dp),dimension(sum(RTI%nlive),sum(RTI%nlive)) :: distance2_matrix
+        real(dp),dimension(settings%nDims,sum(RTI%nlive)) :: position_matrix
         integer,dimension(sum(RTI%nlive)) :: clusters
 
         integer :: num_clusters
         integer :: num_old_clusters
 
         ! number of live points
         integer :: nlive
+        ! number of dimensions
+        integer :: nDims
 
         integer :: i_cluster
 
@@ -295,24 +297,30 @@ function cluster(distance2_matrix) result(cluster_list)
         ! Get the number of old clusters
         num_old_clusters=RTI%ncluster
 
+        ! get the number of dimensions
+        nDims = settings%nDims
+
         i_cluster = 1
         do while(i_cluster<=num_old_clusters)
 
             nlive = RTI%nlive(i_cluster) ! Get the number of live points in a temp variable
 
             if(nlive>2) then 
-                ! Calculate the similarity matrix for this cluster
+                ! get the position matrix for this cluster
                 if(present(sub_dimensions)) then
-                    distance2_matrix(:nlive,:nlive) =&
-                        calculate_similarity_matrix(RTI%live(sub_dimensions,:nlive,i_cluster))
+                    position_matrix(:nDims, :nlive) =&
+                        RTI%live(sub_dimensions,:nlive,i_cluster)
                 else 
-                    distance2_matrix(:nlive,:nlive) =&
-                        calculate_similarity_matrix(RTI%live(settings%h0:settings%h1,:nlive,i_cluster))
+                    position_matrix(:nDims, :nlive) =&
+                        RTI%live(settings%h0:settings%h1,:nlive,i_cluster)
                 end if
 
-                clusters(:nlive) = cluster(distance2_matrix(:nlive,:nlive))
+                clusters(:nlive) = cluster(position_matrix(:nDims, :nlive))
+
+                ! default to KNN clustering
                 if (any(clusters(:nlive)<=0)) then
-                    clusters(:nlive) = NN_clustering(distance2_matrix(:nlive,:nlive))
+                    clusters(:nlive) = NN_clustering(&
+                        calculate_distance2_matrix(position_matrix(:nDims, :nlive)))
                 end if
                 num_clusters = maxval(clusters(:nlive))
 

src/polychord/interfaces.hpp

@@ -47,7 +47,7 @@ void run_polychord(
         double (*c_loglikelihood_ptr)(double*,int,double*,int), 
         void (*c_prior_ptr)(double*,double*,int), 
         void (*c_dumper_ptr)(int,int,int,double*,double*,double*,double,double), 
-        void (*c_cluster_ptr)(double*,int*,int), 
+        void (*c_cluster_ptr)(double*,int*,int,int), 
         Settings s);
 void run_polychord( 
         double (*c_loglikelihood_ptr)(double*,int,double*,int), 
@@ -75,7 +75,7 @@ void run_polychord(
         double (*loglikelihood)(double*,int,double*,int), 
         void (*prior)(double*,double*,int), 
         void (*dumper)(int,int,int,double*,double*,double*,double,double), 
-        void (*c_cluster_ptr)(double*,int*,int), 
+        void (*c_cluster_ptr)(double*,int*,int,int), 
         Settings, MPI_Comm &comm);
 void run_polychord( 
         double (*loglikelihood)(double*,int,double*,int), 
@@ -102,4 +102,4 @@ void run_polychord(
 double default_loglikelihood(double*,int,double*,int); 
 void default_prior(double*,double*,int); 
 void default_dumper(int,int,int,double*,double*,double*,double,double);
-void default_cluster(double*,int*,int);
+void default_cluster(double*,int*,int,int);