swapped indices over

pypolychord/_pypolychord.cpp

@@ -120,7 +120,7 @@ static PyObject *python_cluster = NULL;
 void cluster(double* points, int* cluster_list, int nDims, int nPoints)
 {
     /* create a python version of points */
-    npy_intp shape[] = {nDims,nPoints};            
+    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);

pypolychord/polychord.py

@@ -13,7 +13,7 @@ def default_dumper(live, dead, logweights, logZ, logZerr):
 
 
 def default_cluster(points):
-    return np.zeros(points.shape[1],dtype=int)
+    return np.zeros(points.shape[0],dtype=int)
 
 
 def run_polychord(loglikelihood, nDims, nDerived, settings,
@@ -113,7 +113,7 @@ def run_polychord(loglikelihood, nDims, nDerived, settings,
         Parameters
         ----------
         points: numpy.array
-            positions of points. Shape (nDims, nPoints)
+            positions of points. Shape (nPoints, nDims)
 
         Returns
         -------

run_pypolychord.py

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

src/polychord/c_interface.cpp

@@ -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* points, int* cluster_list, int m, 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,32 +79,32 @@ 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
+    !! Assume that the points can be considered an indexed array of vectors
     !! V = ( v_i : i=1,n )
     !!
-    !! The similarity matrix can be expressed very neatly as
+    !! The distance^2 matrix can be expressed very neatly as
     !! d_ij = (v_i-v_j) . (v_i-v_j)
     !!      = v_i.v_i + v_j.v_j - 2 v_i.v_j
     !!
-    !! The final term can be written as a data_array^T data_array, and the first
+    !! The final term can be written as a points^T points, and the first
     !! two are easy to write. We can therefore calculate this in two lines with
     !! instrisic functions
-    function calculate_distance2_matrix(data_array) result(distance2_matrix)
+    function calculate_distance2_matrix(points) result(distance2_matrix)
 
-        real(dp), intent(in), dimension(:,:) :: data_array
+        real(dp), intent(in), dimension(:,:) :: points
 
-        real(dp), dimension(size(data_array,2),size(data_array,2)) ::distance2_matrix 
+        real(dp), dimension(size(points,1),size(points,1)) ::distance2_matrix 
 
         integer :: i
 
 
         distance2_matrix = spread( &
-            [ ( dot_product(data_array(:,i),data_array(:,i)), i=1,size(data_array,2) ) ], &
-            dim=2,ncopies=size(data_array,2) )
+            [ ( dot_product(points(i,:),points(i,:)), i=1,size(points,1) ) ], &
+            dim=2,ncopies=size(points,1) )
 
-        distance2_matrix = distance2_matrix + transpose(distance2_matrix) - 2d0 * matmul( transpose(data_array),data_array )
+        distance2_matrix = distance2_matrix + transpose(distance2_matrix) - 2d0 * matmul(points,transpose(points))
 
     end function calculate_distance2_matrix
 

src/polychord/clustering.f90

@@ -269,7 +269,7 @@ function do_clustering(settings,RTI,cluster,sub_dimensions)
             function cluster(points) result(cluster_list)
                 import :: dp
                 real(dp), intent(in), dimension(:,:) :: points
-                integer, dimension(size(points,2)) :: cluster_list
+                integer, dimension(size(points,1)) :: cluster_list
             end function
         end interface
 
@@ -308,19 +308,19 @@ function cluster(points) result(cluster_list)
             if(nlive>2) then 
                 ! get the position matrix for this cluster
                 if(present(sub_dimensions)) then
-                    points(:nDims, :nlive) =&
-                        RTI%live(sub_dimensions,:nlive,i_cluster)
+                    points(:nlive, :nDims) =&
+                        transpose(RTI%live(sub_dimensions,:nlive,i_cluster))
                 else 
-                    points(:nDims, :nlive) =&
-                        RTI%live(settings%h0:settings%h1,:nlive,i_cluster)
+                    points(:nlive, :nDims) =&
+                        transpose(RTI%live(settings%h0:settings%h1,:nlive,i_cluster))
                 end if
 
-                clusters(:nlive) = cluster(points(:nDims, :nlive))
+                clusters(:nlive) = cluster(points(:nlive, :nDims))
 
                 ! default to KNN clustering
                 if (any(clusters(:nlive)<=0)) then
                     clusters(:nlive) = NN_clustering(&
-                        calculate_distance2_matrix(points(:nDims, :nlive)))
+                        calculate_distance2_matrix(points(:nlive, :nDims)))
                 end if
                 num_clusters = maxval(clusters(:nlive))
 

src/polychord/interfaces.F90

@@ -50,7 +50,7 @@ end subroutine dumper
             function cluster(points) result(cluster_list)
                 import :: dp
                 real(dp), intent(in), dimension(:,:) :: points
-                integer, dimension(size(points,2)) :: cluster_list
+                integer, dimension(size(points,1)) :: cluster_list
             end function
         end interface
 
@@ -132,7 +132,7 @@ end subroutine dumper
     contains
         function cluster(points) result(cluster_list)
             real(dp), intent(in), dimension(:,:) :: points
-            integer, dimension(size(points,2)) :: cluster_list
+            integer, dimension(size(points,1)) :: cluster_list
             cluster_list = 0
         end function
     end subroutine run_polychord_no_cluster
@@ -390,11 +390,11 @@ subroutine c_dumper(ndead, nlive, npars, live, dead, logweights, logZ, logZerr)
             end subroutine c_dumper
         end interface
         interface
-            subroutine c_cluster(points,cluster_list,m,n) bind(c)
+            subroutine c_cluster(points,cluster_list,nDims,nPoints) bind(c)
                 use iso_c_binding
-                integer(c_int), intent(in), value :: m, n
-                real(c_double), intent(in),  dimension(m,n) :: points
-                integer(c_int), intent(out), dimension(n) :: cluster_list
+                integer(c_int), intent(in), value :: nDims, nPoints
+                real(c_double), intent(in),  dimension(nPoints,nDims) :: points
+                integer(c_int), intent(out), dimension(nPoints) :: cluster_list
             end subroutine c_cluster
         end interface
 
@@ -564,12 +564,12 @@ end subroutine dumper
         function cluster(points) result(cluster_list)
             implicit none
             real(dp), intent(in), dimension(:,:) :: points
-            integer, dimension(size(points,2)) :: cluster_list
+            integer, dimension(size(points,1)) :: cluster_list
 
             integer(c_int) :: c_npoints, c_ndims
             real(c_double),  dimension(size(points,1),size(points,2)) :: c_points
-            integer(c_int), dimension(size(points,2)) :: c_cluster_list
-            c_ndims = size(points,1)
+            integer(c_int), dimension(size(points,1)) :: c_cluster_list
+            c_ndims = size(points,2)
             c_npoints = size(c_cluster_list)
             c_points = points
             call f_cluster_ptr(c_points,c_cluster_list,c_ndims,c_npoints)

src/polychord/nested_sampling.F90

@@ -60,7 +60,7 @@ end subroutine dumper
             function cluster(points) result(cluster_list)
                 import :: dp
                 real(dp), intent(in), dimension(:,:) :: points
-                integer, dimension(size(points,2)) :: cluster_list
+                integer, dimension(size(points,1)) :: cluster_list
             end function
         end interface