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

.github/workflows/CI.yml

@@ -4,7 +4,6 @@ on:
   push:
     branches: [ master ]
   pull_request:
-    branches: [ master ]
   schedule:
     - cron: "0 0 * * *"
 

Makefile_cray

@@ -10,13 +10,13 @@ LDSHARED = $(LD) -shared
 # Archive tool
 AR = ar rv
 
-FFLAGS += -fpp -fpic -qopenmp -dynamic
-CXXFLAGS += -fpic -qopenmp -dynamic
+FFLAGS += -cpp -fpic -fopenmp -dynamic
+CXXFLAGS += -fpic -fopenmp -dynamic
 CFLAGS += -fpic
 #setting for Edison/Cori on NERSC
 ifneq ($(NERSC_HOST),)
-FFLAGS += -axMIC-AVX512,CORE-AVX2
-CXXFLAGS += -axMIC-AVX512,CORE-AVX2
+FFLAGS += -ffree-line-length-none -fallow-argument-mismatch 
+CXXFLAGS += -ffree-line-length-none -fallow-argument-mismatch 
 endif
 
 ifdef IPO

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* points, int* cluster_list, int nDims, int nPoints)
 {
-    /* 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 points */
+    npy_intp shape[] = {nPoints,nDims};            
+    PyObject *array_points = PyArray_SimpleNewFromData(2, shape, NPY_DOUBLE, points);
+    if (array_points ==NULL) throw PythonException();
+    PyArray_CLEARFLAGS(reinterpret_cast<PyArrayObject*>(array_points), NPY_ARRAY_WRITEABLE);
 
     /* create a python version of cluster_list */
-    npy_intp shape1[] = {n};            
+    npy_intp shape1[] = {nPoints};            
     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_points); throw PythonException();}
 
     /* Compute cluster_list from the cluster */
-    PyObject_CallFunctionObjArgs(python_cluster,array_distance2_matrix,array_cluster_list,NULL);
+    PyObject_CallFunctionObjArgs(python_cluster,array_points,array_cluster_list,NULL);
 
     /* Garbage collect */
-    Py_DECREF(array_cluster_list); Py_DECREF(array_distance2_matrix);
+    Py_DECREF(array_cluster_list); Py_DECREF(array_points);
 }
 
 

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(points):
+    return np.full(points.shape[0], -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)
+        points: numpy.array
+            positions of points. Shape (nDims, nPoints)
 
         Returns
         -------
@@ -197,8 +197,11 @@ 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(points, cluster_list):
+        cluster_list[:] = cluster(points) + 1
+
+    settings.grade_dims = [int(d) for d in settings.grade_dims]
+    settings.nlives = {float(logL):int(nlive) for logL, nlive in settings.nlives.items()}
 
     # Run polychord from module library
     _pypolychord.run(wrap_loglikelihood,

run_pypolychord.py

@@ -41,9 +41,9 @@ def dumper(live, dead, logweights, logZ, logZerr):
 
 #| Optional cluster function allow user-defined clustering
 
-def cluster(distance2_matrix):
-    npoints = distance2_matrix.shape[0]
-    clusters = np.ones(npoints, dtype=int)
+def cluster(points):
+    npoints = points.shape[0]
+    clusters = np.full(npoints, -1, dtype=int)
 
     # <do some clustering algorithm to assign clusters>
     # - clusters should be an array of cluster labels for each point

src/polychord/c_interface.cpp

@@ -46,14 +46,14 @@ 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, MPI_Comm& comm)
 #else
 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)
 #endif
 {
@@ -125,7 +125,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)
 {
 	int flag;
@@ -228,5 +228,5 @@ void default_prior(double* cube, double* theta, int nDims)
 
 void default_dumper(int,int,int,double*,double*,double*,double,double) {}
 
-void default_cluster(double* distance2_matrix, int* cluster_list, int n)
-{ for(int i=0;i<n;i++) cluster_list[i] = 0; }
+void default_cluster(double* points, int* cluster_list, int nDims, int nPoints)
+{ for(int i=0;i<nPoints;i++) cluster_list[i] = 0; }

src/polychord/calculate.f90

@@ -79,7 +79,7 @@ function calculate_posterior_point(settings,point,logweight,evidence,volume) res
     end function calculate_posterior_point
 
 
-    !> This function computes the similarity matrix of an array of data.
+    !> This function computes the distance^2 matrix of an array of data.
     !!
     !! Assume that the data_array can be considered an indexed array of vectors
     !! V = ( v_i : i=1,n )
@@ -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

@@ -4,13 +4,13 @@ module KNN_clustering
     implicit none
     contains
 
-    !> This function returns a clustering from a similarity matrix based on
-    !! 'nearest neighbor' clustering.
+    !> This function returns a clustering from a distance^2 matrix based on
+    !! 'nearest neighbour' clustering.
     !!
     !! Points belong to the same cluster if they are in either of each others k
-    !! nearest neighbor sets. 
+    !! nearest neighbour sets. 
     !!
-    !! The algorithm computes the k nearest neihbor sets from the similarity
+    !! The algorithm computes the k nearest neighbor sets from the distance^2
     !! matrix, and then tests
     recursive function NN_clustering(distance2_matrix) result(cluster_list)
         use utils_module, only: relabel
@@ -43,7 +43,7 @@ recursive function NN_clustering(distance2_matrix) result(cluster_list)
         ! 10 degrees of separation are usually fine, we'll expand this if necessary
         k = min(nlive,10)
 
-        ! compute the k nearest neighbors for each point
+        ! compute the k nearest neighbours for each point
         knn(:k,:) = compute_knn(distance2_matrix,k)
 
         ! Set up the cluster list
@@ -83,7 +83,7 @@ recursive function NN_clustering(distance2_matrix) result(cluster_list)
                 ! Get the indices of cluster i_cluster
                 call get_indices_of_cluster(cluster_list,points,i_cluster)
 
-                ! Call this function again on the similarity sub matrix, adding an offset
+                ! Call this function again on the distance^2 sub matrix, adding an offset
                 cluster_list(points) = num_clusters + NN_clustering(distance2_matrix(points,points))
                 num_clusters_new = maxval(cluster_list(points))-num_clusters
 
@@ -114,8 +114,8 @@ function do_clustering_k(knn) result(c)
 
                 ! If they're not in the same cluster already...
                 if(c(i)/=c(j)) then
-                    ! ... check to see if they are within each others k nearest neihbors...
-                    if( neighbors( knn(:,i),knn(:,j) ) ) then
+                    ! ... check to see if they are within each others k nearest neighbours...
+                    if( neighbours( knn(:,i),knn(:,j) ) ) then
 
                         ! If they are then relabel cluster_i and cluster_j to the smaller of the two
                         where(c==c(i).or.c==c(j)) 
@@ -138,10 +138,10 @@ function compute_knn(distance2_matrix,k) result(knn)
         !> The data to compute on
         real(dp), intent(in),dimension(:,:) :: distance2_matrix
 
-        !> The number of nearest neighbors to compute
+        !> The number of nearest neighbours to compute
         integer, intent(in) :: k
 
-        ! The indices of the k nearest neighbors to output
+        ! The indices of the k nearest neighbours to output
         integer, dimension(k,size(distance2_matrix,1)) :: knn
 
         integer :: nPoints,i,j
@@ -155,7 +155,7 @@ function compute_knn(distance2_matrix,k) result(knn)
         knn=0
 
         do i=1,nPoints
-            ! Find the k nearest neighbors for each point
+            ! Find the k nearest neighbours for each point
             distance2s = huge(1d0)
             do j=1,nPoints
                 ! If the distance between i and j is too large to be considered,
@@ -175,22 +175,22 @@ function compute_knn(distance2_matrix,k) result(knn)
     end function compute_knn
 
 
-    ! Return whether they're each others n nearest neighbor list
-    function neighbors(knn1,knn2) result(same_list)
+    ! Return whether they're each others n nearest neighbour list
+    function neighbours(knn1,knn2) result(same_list)
         implicit none
         integer,intent(in), dimension(:) :: knn1
         integer,intent(in), dimension(:) :: knn2 
 
         logical :: same_list
 
-        ! Check to see if they're in each others neighbors lists
+        ! Check to see if they're in each others neighbours lists
         same_list= any(knn1==knn2(1)) .or. any(knn2==knn1(1))
 
-    end function neighbors
+    end function neighbours
 
 
 
-    ! Return the number of matches in the n nearest neighbor list
+    ! Return the number of matches in the n nearest neighbour list
     function matches(knn1,knn2)
         implicit none
         integer,intent(in), dimension(:) :: knn1
@@ -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,26 +266,28 @@ 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(points) 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(:,:) :: points
+                integer, dimension(size(points,2)) :: cluster_list
             end function
         end interface
 
         !> Whether or not a cluster has been found
         logical :: do_clustering
 
 
-        ! Similarity matrix
-        real(dp),dimension(sum(RTI%nlive),sum(RTI%nlive)) :: distance2_matrix
+        ! distance^2 matrix
+        real(dp),dimension(settings%nDims,sum(RTI%nlive)) :: points
         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))
+                    points(: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))
+                    points(:nDims, :nlive) =&
+                        RTI%live(settings%h0:settings%h1,:nlive,i_cluster)
                 end if
 
-                clusters(:nlive) = cluster(distance2_matrix(:nlive,:nlive))
+                clusters(:nlive) = cluster(points(: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(points(:nDims, :nlive)))
                 end if
                 num_clusters = maxval(clusters(:nlive))