From e4622f5b8c3d4b67f04e0fb3baa772a326678778 Mon Sep 17 00:00:00 2001 From: Ormorod Date: Thu, 26 May 2022 19:38:19 +0100 Subject: [PATCH] swapped indices over --- pypolychord/_pypolychord.cpp | 2 +- pypolychord/polychord.py | 4 ++-- run_pypolychord.py | 2 +- src/polychord/c_interface.cpp | 4 ++-- src/polychord/calculate.f90 | 20 ++++++++++---------- src/polychord/clustering.f90 | 14 +++++++------- src/polychord/interfaces.F90 | 18 +++++++++--------- src/polychord/nested_sampling.F90 | 2 +- 8 files changed, 33 insertions(+), 33 deletions(-) diff --git a/pypolychord/_pypolychord.cpp b/pypolychord/_pypolychord.cpp index 1c20a9ca..63d0b43c 100644 --- a/pypolychord/_pypolychord.cpp +++ b/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(array_points), NPY_ARRAY_WRITEABLE); diff --git a/pypolychord/polychord.py b/pypolychord/polychord.py index d916b3f8..7a438814 100644 --- a/pypolychord/polychord.py +++ b/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 ------- diff --git a/run_pypolychord.py b/run_pypolychord.py index c63bff6c..56e05aa8 100755 --- a/run_pypolychord.py +++ b/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) # diff --git a/src/polychord/c_interface.cpp b/src/polychord/c_interface.cpp index 17bf0a01..6cf36d9e 100644 --- a/src/polychord/c_interface.cpp +++ b/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 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 diff --git a/src/polychord/clustering.f90 b/src/polychord/clustering.f90 index 03efdd85..71f81f5c 100644 --- a/src/polychord/clustering.f90 +++ b/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)) diff --git a/src/polychord/interfaces.F90 b/src/polychord/interfaces.F90 index 142aee5b..c2f58487 100644 --- a/src/polychord/interfaces.F90 +++ b/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) diff --git a/src/polychord/nested_sampling.F90 b/src/polychord/nested_sampling.F90 index 8b39dc12..66aef4b9 100644 --- a/src/polychord/nested_sampling.F90 +++ b/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