From f1f52edbd5b8a9593342864dfd1d0f25ac0c8c49 Mon Sep 17 00:00:00 2001 From: Ormorod Date: Thu, 19 May 2022 17:15:31 +0100 Subject: [PATCH 01/39] neighbor spelling (assuming we're sticking to American spelling!) --- src/polychord/clustering.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/polychord/clustering.f90 b/src/polychord/clustering.f90 index 9c49eb8c..4e805812 100644 --- a/src/polychord/clustering.f90 +++ b/src/polychord/clustering.f90 @@ -10,7 +10,7 @@ module KNN_clustering !! Points belong to the same cluster if they are in either of each others k !! nearest neighbor sets. !! - !! The algorithm computes the k nearest neihbor sets from the similarity + !! The algorithm computes the k nearest neighbor sets from the similarity !! matrix, and then tests recursive function NN_clustering(similarity_matrix,num_clusters) result(cluster_list) use utils_module, only: relabel @@ -113,7 +113,7 @@ 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... + ! ... check to see if they are within each others k nearest neighbors... if( neighbors( knn(:,i),knn(:,j) ) ) then ! If they are then relabel cluster_i and cluster_j to the smaller of the two From 1c40b21e6f5043fe664f0552cda258882fa42b90 Mon Sep 17 00:00:00 2001 From: Ormorod Date: Mon, 23 May 2022 17:53:36 +0100 Subject: [PATCH 02/39] custom clustering now takes position matrix instead of distance^2 matrix, to allow more general clustering algorithms --- pypolychord/_pypolychord.cpp | 16 ++++++++-------- pypolychord/polychord.py | 12 ++++++------ run_pypolychord.py | 4 ++-- src/polychord/clustering.f90 | 29 ++++++++++++++++++----------- 4 files changed, 34 insertions(+), 27 deletions(-) diff --git a/pypolychord/_pypolychord.cpp b/pypolychord/_pypolychord.cpp index 03b3794a..eee2b97c 100644 --- a/pypolychord/_pypolychord.cpp +++ b/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 n) { - /* create a python version of distance2_matrix */ + /* create a python version of position_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(array_distance2_matrix), NPY_ARRAY_WRITEABLE); + PyObject *array_position_matrix = PyArray_SimpleNewFromData(2, shape, NPY_DOUBLE, position_matrix); + if (array_position_matrix ==NULL) throw PythonException(); + PyArray_CLEARFLAGS(reinterpret_cast(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); } diff --git a/pypolychord/polychord.py b/pypolychord/polychord.py index 388d5996..9a2b087e 100644 --- a/pypolychord/polychord.py +++ b/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[0],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 (nPoints, nDims) 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, diff --git a/run_pypolychord.py b/run_pypolychord.py index 5aaaf890..1c37e426 100755 --- a/run_pypolychord.py +++ b/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[0] clusters = np.ones(npoints, dtype=int) # diff --git a/src/polychord/clustering.f90 b/src/polychord/clustering.f90 index a2efae04..53eb9d6a 100644 --- a/src/polychord/clustering.f90 +++ b/src/polychord/clustering.f90 @@ -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,7 +278,7 @@ function cluster(distance2_matrix) result(cluster_list) ! Similarity matrix - real(dp),dimension(sum(RTI%nlive),sum(RTI%nlive)) :: distance2_matrix + real(dp),dimension(sum(RTI%nlive),sum(RTI%nlive)) :: position_matrix integer,dimension(sum(RTI%nlive)) :: clusters integer :: num_clusters @@ -286,6 +286,7 @@ function cluster(distance2_matrix) result(cluster_list) ! number of live points integer :: nlive + integer :: nDims integer :: i_cluster @@ -295,24 +296,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(:nlive,:nDims) =& + 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(:nlive,:nDims) =& + RTI%live(settings%h0:settings%h1,:nlive,i_cluster) end if - clusters(:nlive) = cluster(distance2_matrix(:nlive,:nlive)) + clusters(:nlive) = cluster(position_matrix(:nlive,:nDims)) + + ! default to KNN clustering if (any(clusters(:nlive)<=0)) then - clusters(:nlive) = NN_clustering(distance2_matrix(:nlive,:nlive)) + clusters(:nlive) = NN_clustering(& + calculate_similarity_matrix(position_matrix(:nlive, :nDims))) end if num_clusters = maxval(clusters(:nlive)) From 634a4b86eaf101ec48ba391951d159cc61831658 Mon Sep 17 00:00:00 2001 From: Ormorod Date: Tue, 24 May 2022 14:30:57 +0100 Subject: [PATCH 03/39] corrected order of dimensions of position_matrix --- pypolychord/polychord.py | 2 +- run_pypolychord.py | 2 +- src/polychord/clustering.f90 | 11 ++++++----- 3 files changed, 8 insertions(+), 7 deletions(-) diff --git a/pypolychord/polychord.py b/pypolychord/polychord.py index 9a2b087e..10b5c75f 100644 --- a/pypolychord/polychord.py +++ b/pypolychord/polychord.py @@ -113,7 +113,7 @@ def run_polychord(loglikelihood, nDims, nDerived, settings, Parameters ---------- position_matrix: numpy.array - positions of points. Shape (nPoints, nDims) + positions of points. Shape (nDims, nPoints) Returns ------- diff --git a/run_pypolychord.py b/run_pypolychord.py index 1c37e426..8a138068 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(position_matrix): - npoints = position_matrix.shape[0] + npoints = position_matrix.shape[1] clusters = np.ones(npoints, dtype=int) # diff --git a/src/polychord/clustering.f90 b/src/polychord/clustering.f90 index 53eb9d6a..2a70aae2 100644 --- a/src/polychord/clustering.f90 +++ b/src/polychord/clustering.f90 @@ -278,7 +278,7 @@ function cluster(position_matrix) result(cluster_list) ! Similarity matrix - real(dp),dimension(sum(RTI%nlive),sum(RTI%nlive)) :: position_matrix + real(dp),dimension(settings%nDims,sum(RTI%nlive)) :: position_matrix integer,dimension(sum(RTI%nlive)) :: clusters integer :: num_clusters @@ -286,6 +286,7 @@ function cluster(position_matrix) result(cluster_list) ! number of live points integer :: nlive + ! number of dimensions integer :: nDims integer :: i_cluster @@ -307,19 +308,19 @@ function cluster(position_matrix) result(cluster_list) if(nlive>2) then ! get the position matrix for this cluster if(present(sub_dimensions)) then - position_matrix(:nlive,:nDims) =& + position_matrix(:nDims, :nlive) =& RTI%live(sub_dimensions,:nlive,i_cluster) else - position_matrix(:nlive,:nDims) =& + position_matrix(:nDims, :nlive) =& RTI%live(settings%h0:settings%h1,:nlive,i_cluster) end if - clusters(:nlive) = cluster(position_matrix(:nlive,:nDims)) + clusters(:nlive) = cluster(position_matrix(:nDims, :nlive)) ! default to KNN clustering if (any(clusters(:nlive)<=0)) then clusters(:nlive) = NN_clustering(& - calculate_similarity_matrix(position_matrix(:nlive, :nDims))) + calculate_similarity_matrix(position_matrix(:nDims, :nlive))) end if num_clusters = maxval(clusters(:nlive)) From d98cd08b5c31ae278d24d2c0602ad5480d938072 Mon Sep 17 00:00:00 2001 From: Ormorod Date: Tue, 24 May 2022 14:34:51 +0100 Subject: [PATCH 04/39] renamed calculate_similarity_matrix() to calculate_position_matrix() --- src/polychord/calculate.f90 | 4 ++-- src/polychord/clustering.f90 | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/polychord/calculate.f90 b/src/polychord/calculate.f90 index 586e0a62..eca85d54 100644 --- a/src/polychord/calculate.f90 +++ b/src/polychord/calculate.f90 @@ -91,7 +91,7 @@ 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_distance_matrix(data_array) result(similarity_matrix) real(dp), intent(in), dimension(:,:) :: data_array @@ -106,7 +106,7 @@ function calculate_similarity_matrix(data_array) result(similarity_matrix) similarity_matrix = similarity_matrix + transpose(similarity_matrix) - 2d0 * matmul( transpose(data_array),data_array ) - end function calculate_similarity_matrix + end function calculate_distance_matrix diff --git a/src/polychord/clustering.f90 b/src/polychord/clustering.f90 index 2a70aae2..31151fda 100644 --- a/src/polychord/clustering.f90 +++ b/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_distance_matrix use KNN_clustering, only: NN_clustering implicit none @@ -320,7 +320,7 @@ function cluster(position_matrix) result(cluster_list) ! default to KNN clustering if (any(clusters(:nlive)<=0)) then clusters(:nlive) = NN_clustering(& - calculate_similarity_matrix(position_matrix(:nDims, :nlive))) + calculate_distance_matrix(position_matrix(:nDims, :nlive))) end if num_clusters = maxval(clusters(:nlive)) From a04231a27928a46f24fec75f59cef5bc1bd4c176 Mon Sep 17 00:00:00 2001 From: Ormorod Date: Tue, 24 May 2022 14:34:51 +0100 Subject: [PATCH 05/39] renamed calculate_similarity_matrix() to calculate_position2_matrix() --- src/polychord/calculate.f90 | 10 +++++----- src/polychord/clustering.f90 | 4 ++-- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/polychord/calculate.f90 b/src/polychord/calculate.f90 index 586e0a62..941e4011 100644 --- a/src/polychord/calculate.f90 +++ b/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 diff --git a/src/polychord/clustering.f90 b/src/polychord/clustering.f90 index 2a70aae2..2e6a8a13 100644 --- a/src/polychord/clustering.f90 +++ b/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 @@ -320,7 +320,7 @@ function cluster(position_matrix) result(cluster_list) ! default to KNN clustering if (any(clusters(:nlive)<=0)) then clusters(:nlive) = NN_clustering(& - calculate_similarity_matrix(position_matrix(:nDims, :nlive))) + calculate_distance2_matrix(position_matrix(:nDims, :nlive))) end if num_clusters = maxval(clusters(:nlive)) From 956898577e1c349db8fc92365af39545abc7a969 Mon Sep 17 00:00:00 2001 From: Ormorod Date: Wed, 25 May 2022 10:01:32 +0100 Subject: [PATCH 06/39] corrected index in default_cluster --- pypolychord/polychord.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pypolychord/polychord.py b/pypolychord/polychord.py index 10b5c75f..33f7f299 100644 --- a/pypolychord/polychord.py +++ b/pypolychord/polychord.py @@ -13,7 +13,7 @@ def default_dumper(live, dead, logweights, logZ, logZerr): def default_cluster(position_matrix): - return np.zeros(position_matrix.shape[0],dtype=int) + return np.zeros(position_matrix.shape[1],dtype=int) def run_polychord(loglikelihood, nDims, nDerived, settings, From 92b38093d0bf17e79d52d115ed998df5674621cb Mon Sep 17 00:00:00 2001 From: Ormorod Date: Thu, 26 May 2022 15:06:44 +0100 Subject: [PATCH 07/39] added new dimension to cluster input as position_matrix is not square. I believe this is consistent with m = number of dimensions, n = number of clusters --- pypolychord/_pypolychord.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pypolychord/_pypolychord.cpp b/pypolychord/_pypolychord.cpp index eee2b97c..ed709bd0 100644 --- a/pypolychord/_pypolychord.cpp +++ b/pypolychord/_pypolychord.cpp @@ -117,10 +117,10 @@ 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* position_matrix, int* cluster_list, int n) +void cluster(double* position_matrix, int* cluster_list, int m, int n) { /* create a python version of position_matrix */ - npy_intp shape[] = {n,n}; + 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(array_position_matrix), NPY_ARRAY_WRITEABLE); From db02cc29b2d5cf5a10d969ff35b500e0107db879 Mon Sep 17 00:00:00 2001 From: Ormorod Date: Thu, 26 May 2022 15:32:07 +0100 Subject: [PATCH 08/39] changed cluster interfaces to add extra dimension for position_matrix --- src/polychord/interfaces.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/polychord/interfaces.hpp b/src/polychord/interfaces.hpp index 677b4207..17228b8e 100644 --- a/src/polychord/interfaces.hpp +++ b/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); From f5f8c93ff17b35fa17bf490d82dea7498db23d10 Mon Sep 17 00:00:00 2001 From: Will Handley Date: Thu, 26 May 2022 16:02:22 +0100 Subject: [PATCH 09/39] Updated CI to work on pull request to all branches --- .github/workflows/CI.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 4284765e..62bcd96c 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -4,7 +4,6 @@ on: push: branches: [ master ] pull_request: - branches: [ master ] schedule: - cron: "0 0 * * *" From eeb83acb8ba8de65cc8587593e461308d16b7de6 Mon Sep 17 00:00:00 2001 From: Ormorod Date: Thu, 26 May 2022 15:32:07 +0100 Subject: [PATCH 10/39] changed cluster interfaces to add extra dimension for position_matrix --- src/polychord/interfaces.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/polychord/interfaces.hpp b/src/polychord/interfaces.hpp index 677b4207..17228b8e 100644 --- a/src/polychord/interfaces.hpp +++ b/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); From 4d68c9de5a1a6a2bc63d43932273ead511870393 Mon Sep 17 00:00:00 2001 From: Ormorod Date: Thu, 26 May 2022 16:40:12 +0100 Subject: [PATCH 11/39] changed distance2_matrix to no longer be square --- src/polychord/interfaces.F90 | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/src/polychord/interfaces.F90 b/src/polychord/interfaces.F90 index 09dc83b0..f2470652 100644 --- a/src/polychord/interfaces.F90 +++ b/src/polychord/interfaces.F90 @@ -50,7 +50,7 @@ end subroutine dumper function cluster(distance2_matrix) result(cluster_list) import :: dp real(dp), intent(in), dimension(:,:) :: distance2_matrix - integer, dimension(size(distance2_matrix,1)) :: cluster_list + integer, dimension(size(distance2_matrix,2)) :: cluster_list end function end interface @@ -132,7 +132,7 @@ end subroutine dumper contains function cluster(distance2_matrix) result(cluster_list) real(dp), intent(in), dimension(:,:) :: distance2_matrix - integer, dimension(size(distance2_matrix,1)) :: cluster_list + integer, dimension(size(distance2_matrix,2)) :: cluster_list cluster_list = 0 end function end subroutine run_polychord_no_cluster @@ -390,10 +390,10 @@ subroutine c_dumper(ndead, nlive, npars, live, dead, logweights, logZ, logZerr) end subroutine c_dumper end interface interface - subroutine c_cluster(distance2_matrix,cluster_list,n) bind(c) + subroutine c_cluster(distance2_matrix,cluster_list,m,n) bind(c) use iso_c_binding - integer(c_int), intent(in), value :: n - real(c_double), intent(in), dimension(n,n) :: distance2_matrix + integer(c_int), intent(in), value :: m, n + real(c_double), intent(in), dimension(m,n) :: distance2_matrix integer(c_int), intent(out), dimension(n) :: cluster_list end subroutine c_cluster end interface @@ -564,14 +564,15 @@ end subroutine dumper function cluster(distance2_matrix) result(cluster_list) implicit none real(dp), intent(in), dimension(:,:) :: distance2_matrix - integer, dimension(size(distance2_matrix,1)) :: cluster_list + integer, dimension(size(distance2_matrix,2)) :: cluster_list - integer(c_int) :: c_n + integer(c_int) :: c_n, m !m=ndims except I'm scared of the scope real(c_double), dimension(size(distance2_matrix,1),size(distance2_matrix,2)) :: c_distance2_matrix - integer(c_int), dimension(size(distance2_matrix,1)) :: c_cluster_list + integer(c_int), dimension(size(distance2_matrix,2)) :: c_cluster_list + m = size(distance2_matrix,1) c_n = size(c_cluster_list) c_distance2_matrix = distance2_matrix - call f_cluster_ptr(c_distance2_matrix,c_cluster_list,c_n) + call f_cluster_ptr(c_distance2_matrix,c_cluster_list,m,c_n) cluster_list = c_cluster_list end function cluster From b6885cb22fbcbe2c38d92b0ee451d64605f03aee Mon Sep 17 00:00:00 2001 From: Ormorod Date: Thu, 26 May 2022 16:41:35 +0100 Subject: [PATCH 12/39] added extra dimension to cluster argument for non-square distance2_matrix --- src/polychord/c_interface.cpp | 8 ++++---- src/polychord/interfaces.h | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/polychord/c_interface.cpp b/src/polychord/c_interface.cpp index d23b5758..81f5d371 100644 --- a/src/polychord/c_interface.cpp +++ b/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) +void default_cluster(double* distance2_matrix, int* cluster_list, int m, int n) { for(int i=0;i Date: Thu, 26 May 2022 17:29:45 +0100 Subject: [PATCH 13/39] renamed m to nDims, n to nPoints, and position_matrix (and distance2_matrix where still required) to points --- pypolychord/_pypolychord.cpp | 20 ++++++++-------- pypolychord/polychord.py | 10 ++++---- run_pypolychord.py | 4 ++-- src/polychord/c_interface.cpp | 2 +- src/polychord/clustering.f90 | 16 ++++++------- src/polychord/interfaces.F90 | 38 +++++++++++++++---------------- src/polychord/nested_sampling.F90 | 6 ++--- 7 files changed, 48 insertions(+), 48 deletions(-) diff --git a/pypolychord/_pypolychord.cpp b/pypolychord/_pypolychord.cpp index ed709bd0..1c20a9ca 100644 --- a/pypolychord/_pypolychord.cpp +++ b/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* position_matrix, int* cluster_list, int m, int n) +void cluster(double* points, int* cluster_list, int nDims, int nPoints) { - /* 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(array_position_matrix), NPY_ARRAY_WRITEABLE); + /* create a python version of points */ + npy_intp shape[] = {nDims,nPoints}; + 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); /* 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_position_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_position_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_position_matrix); + Py_DECREF(array_cluster_list); Py_DECREF(array_points); } diff --git a/pypolychord/polychord.py b/pypolychord/polychord.py index d24344c9..d916b3f8 100644 --- a/pypolychord/polychord.py +++ b/pypolychord/polychord.py @@ -12,8 +12,8 @@ def default_dumper(live, dead, logweights, logZ, logZerr): pass -def default_cluster(position_matrix): - return np.zeros(position_matrix.shape[1],dtype=int) +def default_cluster(points): + return np.zeros(points.shape[1],dtype=int) def run_polychord(loglikelihood, nDims, nDerived, settings, @@ -112,7 +112,7 @@ def run_polychord(loglikelihood, nDims, nDerived, settings, Parameters ---------- - position_matrix: numpy.array + points: 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(position_matrix, cluster_list): - cluster_list[:] = cluster(position_matrix) + def wrap_cluster(points, cluster_list): + cluster_list[:] = cluster(points) settings.grade_dims = [int(d) for d in settings.grade_dims] settings.nlives = {float(logL):int(nlive) for logL, nlive in settings.nlives.items()} diff --git a/run_pypolychord.py b/run_pypolychord.py index 8a138068..c63bff6c 100755 --- a/run_pypolychord.py +++ b/run_pypolychord.py @@ -41,8 +41,8 @@ def dumper(live, dead, logweights, logZ, logZerr): #| Optional cluster function allow user-defined clustering -def cluster(position_matrix): - npoints = position_matrix.shape[1] +def cluster(points): + npoints = points.shape[1] clusters = np.ones(npoints, dtype=int) # diff --git a/src/polychord/c_interface.cpp b/src/polychord/c_interface.cpp index 81f5d371..17bf0a01 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* distance2_matrix, int* cluster_list, int m, int n) +void default_cluster(double* points, int* cluster_list, int m, int n) { for(int i=0;i2) then ! get the position matrix for this cluster if(present(sub_dimensions)) then - position_matrix(:nDims, :nlive) =& + points(:nDims, :nlive) =& RTI%live(sub_dimensions,:nlive,i_cluster) else - position_matrix(:nDims, :nlive) =& + points(:nDims, :nlive) =& RTI%live(settings%h0:settings%h1,:nlive,i_cluster) end if - clusters(:nlive) = cluster(position_matrix(:nDims, :nlive)) + clusters(:nlive) = cluster(points(:nDims, :nlive)) ! default to KNN clustering if (any(clusters(:nlive)<=0)) then clusters(:nlive) = NN_clustering(& - calculate_distance2_matrix(position_matrix(:nDims, :nlive))) + calculate_distance2_matrix(points(:nDims, :nlive))) end if num_clusters = maxval(clusters(:nlive)) diff --git a/src/polychord/interfaces.F90 b/src/polychord/interfaces.F90 index f2470652..142aee5b 100644 --- a/src/polychord/interfaces.F90 +++ b/src/polychord/interfaces.F90 @@ -47,10 +47,10 @@ subroutine dumper(live, dead, logweights, logZ, logZerr) end subroutine dumper end interface 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,2)) :: cluster_list + real(dp), intent(in), dimension(:,:) :: points + integer, dimension(size(points,2)) :: cluster_list end function end interface @@ -130,9 +130,9 @@ end subroutine dumper #endif call run_polychord(loglikelihood,prior_transform,dumper,cluster,settings,comm) contains - function cluster(distance2_matrix) result(cluster_list) - real(dp), intent(in), dimension(:,:) :: distance2_matrix - integer, dimension(size(distance2_matrix,2)) :: cluster_list + function cluster(points) result(cluster_list) + real(dp), intent(in), dimension(:,:) :: points + integer, dimension(size(points,2)) :: cluster_list cluster_list = 0 end function end subroutine run_polychord_no_cluster @@ -390,10 +390,10 @@ subroutine c_dumper(ndead, nlive, npars, live, dead, logweights, logZ, logZerr) end subroutine c_dumper end interface interface - subroutine c_cluster(distance2_matrix,cluster_list,m,n) bind(c) + subroutine c_cluster(points,cluster_list,m,n) bind(c) use iso_c_binding integer(c_int), intent(in), value :: m, n - real(c_double), intent(in), dimension(m,n) :: distance2_matrix + real(c_double), intent(in), dimension(m,n) :: points integer(c_int), intent(out), dimension(n) :: cluster_list end subroutine c_cluster end interface @@ -561,18 +561,18 @@ subroutine dumper(live, dead, logweights, logZ, logZerr) call f_dumper_ptr(c_ndead, c_nlive, c_npars, c_live, c_dead, c_logweights, c_logZ, c_logZerr) end subroutine dumper - function cluster(distance2_matrix) result(cluster_list) + function cluster(points) result(cluster_list) implicit none - real(dp), intent(in), dimension(:,:) :: distance2_matrix - integer, dimension(size(distance2_matrix,2)) :: cluster_list - - integer(c_int) :: c_n, m !m=ndims except I'm scared of the scope - real(c_double), dimension(size(distance2_matrix,1),size(distance2_matrix,2)) :: c_distance2_matrix - integer(c_int), dimension(size(distance2_matrix,2)) :: c_cluster_list - m = size(distance2_matrix,1) - c_n = size(c_cluster_list) - c_distance2_matrix = distance2_matrix - call f_cluster_ptr(c_distance2_matrix,c_cluster_list,m,c_n) + real(dp), intent(in), dimension(:,:) :: points + integer, dimension(size(points,2)) :: 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) + c_npoints = size(c_cluster_list) + c_points = points + call f_cluster_ptr(c_points,c_cluster_list,c_ndims,c_npoints) cluster_list = c_cluster_list end function cluster diff --git a/src/polychord/nested_sampling.F90 b/src/polychord/nested_sampling.F90 index b903dbe8..66aef4b9 100644 --- a/src/polychord/nested_sampling.F90 +++ b/src/polychord/nested_sampling.F90 @@ -57,10 +57,10 @@ subroutine dumper(live, dead, logweights, logZ, logZerr) end subroutine dumper end interface 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,1)) :: cluster_list end function end interface From fa5c42d55e2860c22be257aa5082d16e1ab44691 Mon Sep 17 00:00:00 2001 From: Ormorod Date: Thu, 26 May 2022 17:33:20 +0100 Subject: [PATCH 14/39] changed neighbour to ~~correct~~ UK spelling --- src/polychord/clustering.f90 | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/polychord/clustering.f90 b/src/polychord/clustering.f90 index 8c7a00e7..c6902621 100644 --- a/src/polychord/clustering.f90 +++ b/src/polychord/clustering.f90 @@ -5,10 +5,10 @@ module KNN_clustering contains !> This function returns a clustering from a similarity matrix based on - !! 'nearest neighbor' clustering. + !! '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 !! matrix, and then tests @@ -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 @@ -115,7 +115,7 @@ 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 + 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 From 51848350adee7b5b5c29c56297f64ed86cf505c7 Mon Sep 17 00:00:00 2001 From: Ormorod Date: Thu, 26 May 2022 18:23:49 +0100 Subject: [PATCH 15/39] forgot to change the interface in clustering itself --- src/polychord/clustering.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/polychord/clustering.f90 b/src/polychord/clustering.f90 index c6902621..03efdd85 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,1)) :: cluster_list + integer, dimension(size(points,2)) :: cluster_list end function end interface From 3ed846cf5f4e564163e0dafe5cfa16f2dc9540a6 Mon Sep 17 00:00:00 2001 From: Ormorod Date: Thu, 26 May 2022 19:33:15 +0100 Subject: [PATCH 16/39] another clustering interface using the wrong dimension of points --- src/polychord/nested_sampling.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/polychord/nested_sampling.F90 b/src/polychord/nested_sampling.F90 index 66aef4b9..8b39dc12 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,1)) :: cluster_list + integer, dimension(size(points,2)) :: cluster_list end function end interface From e4622f5b8c3d4b67f04e0fb3baa772a326678778 Mon Sep 17 00:00:00 2001 From: Ormorod Date: Thu, 26 May 2022 19:38:19 +0100 Subject: [PATCH 17/39] 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 From 23d8db874070226d2e1f89b6fd6a320ca68b63dd Mon Sep 17 00:00:00 2001 From: Ormorod Date: Thu, 26 May 2022 23:38:57 +0100 Subject: [PATCH 18/39] Revert "swapped indices over" This reverts commit e4622f5b8c3d4b67f04e0fb3baa772a326678778. --- 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 63d0b43c..1c20a9ca 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[] = {nPoints,nDims}; + npy_intp shape[] = {nDims,nPoints}; 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 7a438814..d916b3f8 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[0],dtype=int) + return np.zeros(points.shape[1],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 (nPoints, nDims) + positions of points. Shape (nDims, nPoints) Returns ------- diff --git a/run_pypolychord.py b/run_pypolychord.py index 56e05aa8..c63bff6c 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[0] + npoints = points.shape[1] clusters = np.ones(npoints, dtype=int) # diff --git a/src/polychord/c_interface.cpp b/src/polychord/c_interface.cpp index 6cf36d9e..17bf0a01 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 nDims, int nPoints) -{ for(int i=0;i This function computes the distance^2 matrix of an array of data. + !> This function computes the similarity matrix of an array of data. !! - !! Assume that the points can be considered an indexed array of vectors + !! Assume that the data_array can be considered an indexed array of vectors !! V = ( v_i : i=1,n ) !! - !! The distance^2 matrix can be expressed very neatly as + !! The similarity 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 points^T points, and the first + !! 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_distance2_matrix(points) result(distance2_matrix) + function calculate_distance2_matrix(data_array) result(distance2_matrix) - real(dp), intent(in), dimension(:,:) :: points + real(dp), intent(in), dimension(:,:) :: data_array - real(dp), dimension(size(points,1),size(points,1)) ::distance2_matrix + real(dp), dimension(size(data_array,2),size(data_array,2)) ::distance2_matrix integer :: i distance2_matrix = spread( & - [ ( dot_product(points(i,:),points(i,:)), i=1,size(points,1) ) ], & - dim=2,ncopies=size(points,1) ) + [ ( dot_product(data_array(:,i),data_array(:,i)), i=1,size(data_array,2) ) ], & + dim=2,ncopies=size(data_array,2) ) - distance2_matrix = distance2_matrix + transpose(distance2_matrix) - 2d0 * matmul(points,transpose(points)) + distance2_matrix = distance2_matrix + transpose(distance2_matrix) - 2d0 * matmul( transpose(data_array),data_array ) end function calculate_distance2_matrix diff --git a/src/polychord/clustering.f90 b/src/polychord/clustering.f90 index fc79007e..7b1b3d16 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,1)) :: cluster_list + integer, dimension(size(points,2)) :: 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(:nlive, :nDims) =& - transpose(RTI%live(sub_dimensions,:nlive,i_cluster)) + points(:nDims, :nlive) =& + RTI%live(sub_dimensions,:nlive,i_cluster) else - points(:nlive, :nDims) =& - transpose(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(points(:nlive, :nDims)) + clusters(:nlive) = cluster(points(:nDims, :nlive)) ! default to KNN clustering if (any(clusters(:nlive)<=0)) then clusters(:nlive) = NN_clustering(& - calculate_distance2_matrix(points(:nlive, :nDims))) + calculate_distance2_matrix(points(:nDims, :nlive))) end if num_clusters = maxval(clusters(:nlive)) diff --git a/src/polychord/interfaces.F90 b/src/polychord/interfaces.F90 index c2f58487..142aee5b 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,1)) :: cluster_list + integer, dimension(size(points,2)) :: 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,1)) :: cluster_list + integer, dimension(size(points,2)) :: 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,nDims,nPoints) bind(c) + subroutine c_cluster(points,cluster_list,m,n) bind(c) use iso_c_binding - 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 + 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 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,1)) :: cluster_list + integer, dimension(size(points,2)) :: 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,1)) :: c_cluster_list - c_ndims = size(points,2) + integer(c_int), dimension(size(points,2)) :: c_cluster_list + c_ndims = size(points,1) 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 66aef4b9..8b39dc12 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,1)) :: cluster_list + integer, dimension(size(points,2)) :: cluster_list end function end interface From 4a9bf0ac478ec211aae3e93468842b78288450ec Mon Sep 17 00:00:00 2001 From: Ormorod Date: Thu, 26 May 2022 23:55:10 +0100 Subject: [PATCH 19/39] now uses c indexing in c and python, and fortran indexing in fortran --- pypolychord/_pypolychord.cpp | 2 +- pypolychord/polychord.py | 2 +- run_pypolychord.py | 2 +- src/polychord/c_interface.cpp | 4 ++-- src/polychord/calculate.f90 | 2 +- src/polychord/interfaces.F90 | 8 ++++---- 6 files changed, 10 insertions(+), 10 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..0b8cd179 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, 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 !! V = ( v_i : i=1,n ) diff --git a/src/polychord/interfaces.F90 b/src/polychord/interfaces.F90 index 142aee5b..093f359b 100644 --- a/src/polychord/interfaces.F90 +++ b/src/polychord/interfaces.F90 @@ -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(nDims,nPoints) :: points + integer(c_int), intent(out), dimension(nPoints) :: cluster_list end subroutine c_cluster end interface From f68dcb3c553959f4ab82fc3ae3a5b622fb96a65d Mon Sep 17 00:00:00 2001 From: Ormorod Date: Mon, 13 Jun 2022 16:38:07 +0100 Subject: [PATCH 20/39] similarity matrix renamed to distance^2 in comments --- src/polychord/clustering.f90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/polychord/clustering.f90 b/src/polychord/clustering.f90 index 7b1b3d16..c38ba71c 100644 --- a/src/polychord/clustering.f90 +++ b/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 + !> 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 neighbour sets. !! - !! The algorithm computes the k nearest neighbor 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 @@ -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 @@ -277,7 +277,7 @@ function cluster(points) result(cluster_list) logical :: do_clustering - ! Similarity matrix + ! distance^2 matrix real(dp),dimension(settings%nDims,sum(RTI%nlive)) :: points integer,dimension(sum(RTI%nlive)) :: clusters From da786338aef1f7914c447abe73b285bcb2ca2d69 Mon Sep 17 00:00:00 2001 From: Ormorod Date: Thu, 30 Jun 2022 11:05:21 +0100 Subject: [PATCH 21/39] wrap_cluster() adds one to the cluster list, as Python algos number clusters from 0, but fortran wants them numbered from 1 --- pypolychord/polychord.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pypolychord/polychord.py b/pypolychord/polychord.py index 0b8cd179..eea20fba 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[0],dtype=int) + return np.full(points.shape[0], -1,dtype=int) def run_polychord(loglikelihood, nDims, nDerived, settings, @@ -198,7 +198,7 @@ def wrap_prior(cube, theta): theta[:] = prior(cube) def wrap_cluster(points, cluster_list): - cluster_list[:] = cluster(points) + 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()} From be3467a6a98e6e5453614c4cf67bf95c389f8d17 Mon Sep 17 00:00:00 2001 From: Ormorod Date: Thu, 30 Jun 2022 12:05:56 +0100 Subject: [PATCH 22/39] example cluster now returns -1 --- run_pypolychord.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/run_pypolychord.py b/run_pypolychord.py index 56e05aa8..cefdee8d 100755 --- a/run_pypolychord.py +++ b/run_pypolychord.py @@ -43,7 +43,7 @@ def dumper(live, dead, logweights, logZ, logZerr): def cluster(points): npoints = points.shape[0] - clusters = np.ones(npoints, dtype=int) + clusters = np.full(npoints, -1, dtype=int) # # - clusters should be an array of cluster labels for each point From 14efda18a8fce6cdb55476f681c4835225b8b731 Mon Sep 17 00:00:00 2001 From: Ormorod Date: Fri, 29 Jul 2022 14:54:55 +0100 Subject: [PATCH 23/39] example to test custom clustering --- run_pypolychord_custom_clustering.py | 155 +++++++++++++++++++++++++++ 1 file changed, 155 insertions(+) create mode 100644 run_pypolychord_custom_clustering.py diff --git a/run_pypolychord_custom_clustering.py b/run_pypolychord_custom_clustering.py new file mode 100644 index 00000000..9bb149ad --- /dev/null +++ b/run_pypolychord_custom_clustering.py @@ -0,0 +1,155 @@ + +from numpy import pi, log, sqrt +import numpy as np +import pypolychord +from pypolychord.settings import PolyChordSettings +from pypolychord.priors import UniformPrior +from scipy.special import logsumexp +try: + from mpi4py import MPI +except ImportError: + pass + +from sklearn.neighbors import NearestNeighbors + + +#| Define a four-dimensional spherical gaussian likelihood, +#| width sigma=0.1, centered on the 0 with one derived parameter. +#| The derived parameter is the squared radius + +nDims = 4 +nDerived = 1 +sigma = 0.1 + +centres = np.array([[-0.5] * nDims, [0.5] * nDims]) +# weights = np.array([0.5, 0.5]) +weights = np.array([2/3, 1/3]) + + +def likelihood(theta): + """Double Gaussian likelihood.""" + nDims = len(theta) + logL = -np.log(2 * np.pi * sigma * sigma) * nDims / 2 + logL += logsumexp(-np.sum((theta - centres) ** 2, axis=-1) / 2 / sigma / sigma, b = weights) + return logL, [np.sum(theta**2)] + + +#| Define a box uniform prior from -1 to 1 + +def prior(hypercube): + """ Uniform prior from [-1,1]^D. """ + return UniformPrior(-1, 1)(hypercube) + + +#| Optional cluster function allow user-defined clustering +## KNN clustering algorithm translated from clustering.F90 + +def mutual_nearest_neighbors(knn_array, i, ii): + return (i in knn_array[ii]) or (ii in knn_array[i]) + + +def relabel(labels): + k = max(labels) + appearance_order = {} + num_found = 0 + + for label in labels: + if label not in appearance_order: + appearance_order[label] = num_found + num_found += 1 + + for i, old_label in enumerate(labels): + labels[i] = appearance_order[old_label] + return labels + + +def do_clustering(knn_array): + num_points = knn_array.shape[0] + labels = np.arange(num_points) + for iii in range(num_points): + for ii in range(iii + 1, num_points): + if labels[ii] != labels[iii]: + if mutual_nearest_neighbors(knn_array, ii, iii): + for i in range(num_points): + if labels[i] == labels[ii] or labels[i] == labels[iii]: + labels[i] = min([labels[ii], labels[iii]]) + return relabel(labels) + + +def knn_cluster(position_matrix): + """slight concern if two points are the same because sklearn""" + npoints = position_matrix.shape[0] + k = min(npoints, 10) + nn = NearestNeighbors(n_neighbors=k).fit(position_matrix) + knn_array = nn.kneighbors(position_matrix, return_distance=False) + labels = np.arange(npoints) + num_clusters = npoints + + labels_old = labels + num_clusters_old = num_clusters + + for n in np.arange(2, k): + labels = do_clustering(knn_array[:, :n]) + num_clusters = max(labels) + 1 + + if num_clusters <= 0: + raise ValueError("somehow got <= 0 clusters") + elif 1 == num_clusters: + return labels + elif np.all(labels == labels_old): + break + elif k == n - 1: + k = min(k * 2, npoints) + nn = NearestNeighbors(n_neighbors=k).fit(position_matrix) + knn_array = nn.kneighbors(position_matrix, return_distance=False) + + labels_old = labels + num_clusters_old = num_clusters + + if num_clusters > 1: + i_cluster = 0 + while i_cluster < num_clusters: + cluster = position_matrix[labels == i_cluster] + labels[labels == i_cluster] = knn_cluster(cluster) + num_clusters + labels = relabel(labels) + if num_clusters - 1 == max(labels): + i_cluster += 1 + return labels + + +#| Initialise the settings + +settings = PolyChordSettings(nDims, nDerived) +settings.file_root = "custom_clustering" +settings.nlive = 200 +settings.do_clustering = True +settings.read_resume = False + +#| Run PolyChord + +output = pypolychord.run_polychord(likelihood, nDims, nDerived, settings, prior, cluster=knn_cluster) + +#| Create a paramnames file + +paramnames = [('p%i' % i, r'\theta_%i' % i) for i in range(nDims)] +paramnames += [('r*', 'r')] +output.make_paramnames_files(paramnames) + +#| Make an anesthetic plot (could also use getdist) +try: + from anesthetic import NestedSamples + samples = NestedSamples(root= settings.base_dir + '/' + settings.file_root) + fig, axes = samples.plot_2d(['p0','p1','p2','p3','r']) + fig.savefig('posterior.pdf') + +except ImportError: + try: + import getdist.plots + posterior = output.posterior + g = getdist.plots.getSubplotPlotter() + g.triangle_plot(posterior, filled=True) + g.export('custom_clustering.pdf') + except ImportError: + print("Install matplotlib and getdist for plotting examples") + + print("Install anesthetic or getdist for for plotting examples") From 3029e60de19a3d475fc87700777c54aee09533cf Mon Sep 17 00:00:00 2001 From: Ormorod Date: Fri, 29 Jul 2022 15:27:59 +0100 Subject: [PATCH 24/39] equal weights of each Gaussian peak, corrected comments --- run_pypolychord_custom_clustering.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/run_pypolychord_custom_clustering.py b/run_pypolychord_custom_clustering.py index 9bb149ad..b708e593 100644 --- a/run_pypolychord_custom_clustering.py +++ b/run_pypolychord_custom_clustering.py @@ -13,21 +13,20 @@ from sklearn.neighbors import NearestNeighbors -#| Define a four-dimensional spherical gaussian likelihood, -#| width sigma=0.1, centered on the 0 with one derived parameter. -#| The derived parameter is the squared radius +#| Define a four-dimensional twin-spherical gaussian likelihood, +#| width sigma=0.1, centered on ±[0.5]^nDims with one derived parameter. +#| The derived parameter is the squared radius from [0]^nDims nDims = 4 nDerived = 1 sigma = 0.1 centres = np.array([[-0.5] * nDims, [0.5] * nDims]) -# weights = np.array([0.5, 0.5]) -weights = np.array([2/3, 1/3]) +weights = np.array([0.5, 0.5]) def likelihood(theta): - """Double Gaussian likelihood.""" + """Twin Gaussian likelihood.""" nDims = len(theta) logL = -np.log(2 * np.pi * sigma * sigma) * nDims / 2 logL += logsumexp(-np.sum((theta - centres) ** 2, axis=-1) / 2 / sigma / sigma, b = weights) From 88572e33d3284189d2b35c18547d1c28eeaac996 Mon Sep 17 00:00:00 2001 From: Ormorod Date: Fri, 29 Jul 2022 15:33:51 +0100 Subject: [PATCH 25/39] remove mutual_nearest_neighbours function --- run_pypolychord_custom_clustering.py | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/run_pypolychord_custom_clustering.py b/run_pypolychord_custom_clustering.py index b708e593..672aa550 100644 --- a/run_pypolychord_custom_clustering.py +++ b/run_pypolychord_custom_clustering.py @@ -43,10 +43,6 @@ def prior(hypercube): #| Optional cluster function allow user-defined clustering ## KNN clustering algorithm translated from clustering.F90 -def mutual_nearest_neighbors(knn_array, i, ii): - return (i in knn_array[ii]) or (ii in knn_array[i]) - - def relabel(labels): k = max(labels) appearance_order = {} @@ -68,7 +64,7 @@ def do_clustering(knn_array): for iii in range(num_points): for ii in range(iii + 1, num_points): if labels[ii] != labels[iii]: - if mutual_nearest_neighbors(knn_array, ii, iii): + if (ii in knn_array[iii]) or (iii in knn_array[ii]): for i in range(num_points): if labels[i] == labels[ii] or labels[i] == labels[iii]: labels[i] = min([labels[ii], labels[iii]]) From b49fb80a8e4e52f4d261db00dd96c5b3bc0b7c34 Mon Sep 17 00:00:00 2001 From: Ormorod Date: Fri, 29 Jul 2022 16:27:14 +0100 Subject: [PATCH 26/39] replace sklearn dependence with compute_knn --- run_pypolychord_custom_clustering.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/run_pypolychord_custom_clustering.py b/run_pypolychord_custom_clustering.py index 672aa550..ab5c44d9 100644 --- a/run_pypolychord_custom_clustering.py +++ b/run_pypolychord_custom_clustering.py @@ -4,13 +4,13 @@ import pypolychord from pypolychord.settings import PolyChordSettings from pypolychord.priors import UniformPrior +from scipy.spatial import distance_matrix from scipy.special import logsumexp try: from mpi4py import MPI except ImportError: pass -from sklearn.neighbors import NearestNeighbors #| Define a four-dimensional twin-spherical gaussian likelihood, @@ -43,6 +43,10 @@ def prior(hypercube): #| Optional cluster function allow user-defined clustering ## KNN clustering algorithm translated from clustering.F90 +def compute_knn(position_matrix, k): + return np.argsort(distance_matrix(position_matrix, position_matrix))[:, :k] + + def relabel(labels): k = max(labels) appearance_order = {} @@ -75,8 +79,7 @@ def knn_cluster(position_matrix): """slight concern if two points are the same because sklearn""" npoints = position_matrix.shape[0] k = min(npoints, 10) - nn = NearestNeighbors(n_neighbors=k).fit(position_matrix) - knn_array = nn.kneighbors(position_matrix, return_distance=False) + knn_array = compute_knn(position_matrix, k) labels = np.arange(npoints) num_clusters = npoints @@ -95,8 +98,7 @@ def knn_cluster(position_matrix): break elif k == n - 1: k = min(k * 2, npoints) - nn = NearestNeighbors(n_neighbors=k).fit(position_matrix) - knn_array = nn.kneighbors(position_matrix, return_distance=False) + knn_array = compute_knn(position_matrix, k) labels_old = labels num_clusters_old = num_clusters From f8dc320090c18b306c3fb6324ffbef59222d4ec4 Mon Sep 17 00:00:00 2001 From: Ormorod Date: Fri, 29 Jul 2022 16:36:41 +0100 Subject: [PATCH 27/39] compute_knn() computed the full nn ordering, then returned the first k, so instead just return all of them (compute_nn()) and don't bother recalculating. --- run_pypolychord_custom_clustering.py | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/run_pypolychord_custom_clustering.py b/run_pypolychord_custom_clustering.py index ab5c44d9..8804e450 100644 --- a/run_pypolychord_custom_clustering.py +++ b/run_pypolychord_custom_clustering.py @@ -43,8 +43,8 @@ def prior(hypercube): #| Optional cluster function allow user-defined clustering ## KNN clustering algorithm translated from clustering.F90 -def compute_knn(position_matrix, k): - return np.argsort(distance_matrix(position_matrix, position_matrix))[:, :k] +def compute_nn(position_matrix): + return np.argsort(distance_matrix(position_matrix, position_matrix)) def relabel(labels): @@ -79,15 +79,14 @@ def knn_cluster(position_matrix): """slight concern if two points are the same because sklearn""" npoints = position_matrix.shape[0] k = min(npoints, 10) - knn_array = compute_knn(position_matrix, k) + nn_array = compute_nn(position_matrix) labels = np.arange(npoints) num_clusters = npoints labels_old = labels - num_clusters_old = num_clusters for n in np.arange(2, k): - labels = do_clustering(knn_array[:, :n]) + labels = do_clustering(nn_array[:, :n]) num_clusters = max(labels) + 1 if num_clusters <= 0: @@ -98,10 +97,8 @@ def knn_cluster(position_matrix): break elif k == n - 1: k = min(k * 2, npoints) - knn_array = compute_knn(position_matrix, k) labels_old = labels - num_clusters_old = num_clusters if num_clusters > 1: i_cluster = 0 From 4d950776e4ecd91a45128aea96936f55a9e0fc7d Mon Sep 17 00:00:00 2001 From: Ormorod Date: Fri, 29 Jul 2022 16:55:25 +0100 Subject: [PATCH 28/39] realised the for loop in the copied knn clustering is actually a while loop in disguise --- run_pypolychord_custom_clustering.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/run_pypolychord_custom_clustering.py b/run_pypolychord_custom_clustering.py index 8804e450..f3b52413 100644 --- a/run_pypolychord_custom_clustering.py +++ b/run_pypolychord_custom_clustering.py @@ -85,7 +85,8 @@ def knn_cluster(position_matrix): labels_old = labels - for n in np.arange(2, k): + n = 2 + while n <= k: labels = do_clustering(nn_array[:, :n]) num_clusters = max(labels) + 1 @@ -95,10 +96,11 @@ def knn_cluster(position_matrix): return labels elif np.all(labels == labels_old): break - elif k == n - 1: + elif k == n: k = min(k * 2, npoints) labels_old = labels + n += 1 if num_clusters > 1: i_cluster = 0 From 2d99c895d96eed1ad2b2b1847f1bf56a8791c180 Mon Sep 17 00:00:00 2001 From: Ormorod Date: Fri, 29 Jul 2022 16:58:15 +0100 Subject: [PATCH 29/39] Removed unused num_clusters_old from native clustering --- src/polychord/clustering.f90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/polychord/clustering.f90 b/src/polychord/clustering.f90 index c38ba71c..fdfd4ffc 100644 --- a/src/polychord/clustering.f90 +++ b/src/polychord/clustering.f90 @@ -48,7 +48,6 @@ recursive function NN_clustering(distance2_matrix) result(cluster_list) ! Set up the cluster list cluster_list_old = [( i_point,i_point=1,nlive )] - num_clusters_old = nlive do n=2,k @@ -71,7 +70,6 @@ recursive function NN_clustering(distance2_matrix) result(cluster_list) ! Save the old cluster list for later. cluster_list_old = cluster_list - num_clusters_old = num_clusters end do From 1a561de55dc358314b04bfac339e4b758870366a Mon Sep 17 00:00:00 2001 From: Ormorod Date: Fri, 29 Jul 2022 17:02:44 +0100 Subject: [PATCH 30/39] comment reminds that clustering algorithm should be 0-indexed --- run_pypolychord.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/run_pypolychord.py b/run_pypolychord.py index cefdee8d..e03c818b 100755 --- a/run_pypolychord.py +++ b/run_pypolychord.py @@ -48,7 +48,8 @@ def cluster(points): # # - clusters should be an array of cluster labels for each point # - each cluster should have at least one point - # - thus max(clusters) should be the number of clusters + # - thus max(clusters)+1 should be the number of clusters + # - i.e. clusters are 0-indexed # - work with the above numpy integer array return clusters From f67543d0c68a52ccda8523ec5159637518e13d1e Mon Sep 17 00:00:00 2001 From: Ormorod Date: Fri, 29 Jul 2022 17:07:31 +0100 Subject: [PATCH 31/39] updated run_pypolychord.ipynb with py2nb --- run_pypolychord.ipynb | 158 ++++++++++++++++++------------------------ 1 file changed, 66 insertions(+), 92 deletions(-) diff --git a/run_pypolychord.ipynb b/run_pypolychord.ipynb index 79d69f78..31021fb6 100644 --- a/run_pypolychord.ipynb +++ b/run_pypolychord.ipynb @@ -2,18 +2,25 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": null, + "id": "07ce1c6b", "metadata": {}, "outputs": [], "source": [ "from numpy import pi, log, sqrt\n", + "import numpy as np\n", "import pypolychord\n", "from pypolychord.settings import PolyChordSettings\n", - "from pypolychord.priors import UniformPrior" + "from pypolychord.priors import UniformPrior\n", + "try:\n", + " from mpi4py import MPI\n", + "except ImportError:\n", + " pass" ] }, { "cell_type": "markdown", + "id": "ddf5dfd5", "metadata": {}, "source": [ "Define a four-dimensional spherical gaussian likelihood,\n", @@ -23,7 +30,8 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, + "id": "071a0561", "metadata": {}, "outputs": [], "source": [ @@ -44,6 +52,7 @@ }, { "cell_type": "markdown", + "id": "186e760a", "metadata": {}, "source": [ "Define a box uniform prior from -1 to 1" @@ -51,7 +60,8 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": null, + "id": "c6154b8d", "metadata": {}, "outputs": [], "source": [ @@ -62,6 +72,7 @@ }, { "cell_type": "markdown", + "id": "54a00017", "metadata": {}, "source": [ "Optional dumper function giving run-time read access to\n", @@ -70,7 +81,8 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": null, + "id": "9a7c7238", "metadata": {}, "outputs": [], "source": [ @@ -80,6 +92,36 @@ }, { "cell_type": "markdown", + "id": "a93f0471", + "metadata": {}, + "source": [ + "Optional cluster function allow user-defined clustering" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b1b4eaaf", + "metadata": {}, + "outputs": [], + "source": [ + "def cluster(points):\n", + " npoints = points.shape[0]\n", + " clusters = np.full(npoints, -1, dtype=int)\n", + "\n", + " # \n", + " # - clusters should be an array of cluster labels for each point\n", + " # - each cluster should have at least one point\n", + " # - thus max(clusters)+1 should be the number of clusters\n", + " # - i.e. clusters are 0-indexed\n", + " # - work with the above numpy integer array\n", + "\n", + " return clusters" + ] + }, + { + "cell_type": "markdown", + "id": "7ccf5541", "metadata": {}, "source": [ "Initialise the settings" @@ -87,7 +129,8 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": null, + "id": "caf7fb7c", "metadata": {}, "outputs": [], "source": [ @@ -100,6 +143,7 @@ }, { "cell_type": "markdown", + "id": "8ca5896d", "metadata": {}, "source": [ "Run PolyChord" @@ -107,52 +151,17 @@ }, { "cell_type": "code", - "execution_count": 6, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Last dead point: [ -0.12195646 0.37405823 -0.74549651 -0.57059283 1.03613417\n", - " -107.50194883 -46.27212213]\n", - "Last dead point: [ -0.52163002 -0.06336759 -0.43192652 0.38364359 0.60985625\n", - " -41.35375964 -24.95822646]\n", - "Last dead point: [ -0.37403331 0.03735235 -0.32821553 -0.31645229 0.3491636\n", - " -19.55187971 -11.92359365]\n", - "Last dead point: [-0.11967879 0.22564894 0.34498611 0.12669236 0.20030682 -7.93394934\n", - " -4.480755 ]\n", - "Last dead point: [ 0.15947558 -0.17120628 0.21292861 -0.12136559 0.11481225 -0.97442395\n", - " -0.20602622]\n", - "Last dead point: [-0.05449124 -0.24856098 0.053324 -0.05999985 0.07119528 1.95448519\n", - " 1.97482204]\n", - "Last dead point: [ 0.11354188 0.04542128 0.04324486 0.15991583 0.04239804 -0.65204203\n", - " 3.41468411]\n", - "Last dead point: [0.03245751 0.08394396 0.12165265 0.04758947 0.0251642 4.03292303\n", - " 4.2763761 ]\n", - "Last dead point: [ 0.05227755 -0.01961902 -0.06877328 -0.08374121 0.0148602 4.42213521\n", - " 4.79157616]\n", - "Last dead point: [7.75444943e-02 3.43045989e-02 3.54071198e-03 4.29640866e-02\n", - " 9.04840349e-03 4.89388098e+00 5.08216606e+00]\n", - "Last dead point: [0.01163431 0.01282317 0.04719472 0.05280419 0.00531541 5.02547622\n", - " 5.26881552]\n", - "Last dead point: [-4.50972426e-02 2.65906632e-02 -1.05166059e-02 -2.14135125e-02\n", - " 3.30996217e-03 5.35102039e+00 5.36908813e+00]\n", - "Last dead point: [-3.64909631e-03 -1.36859157e-02 6.78829623e-03 -4.23035741e-02\n", - " 2.03629354e-03 5.11415470e+00 5.43277156e+00]\n", - "Last dead point: [-8.32623304e-03 -7.46293828e-03 -2.53157193e-02 1.98799482e-02\n", - " 1.16111959e-03 5.45981438e+00 5.47653026e+00]\n", - "Last dead point: [ 4.69883258e-03 1.81089351e-03 -4.12460698e-03 -2.22366928e-03\n", - " 4.73154508e-05 5.48264395e+00 5.53222047e+00]\n" - ] - } - ], - "source": [ - "output = pypolychord.run_polychord(likelihood, nDims, nDerived, settings, prior, dumper)" + "execution_count": null, + "id": "fcab6df8", + "metadata": {}, + "outputs": [], + "source": [ + "output = pypolychord.run_polychord(likelihood, nDims, nDerived, settings, prior, dumper, cluster)" ] }, { "cell_type": "markdown", + "id": "3b6b8b82", "metadata": {}, "source": [ "Create a paramnames file" @@ -160,7 +169,8 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": null, + "id": "0b20d84f", "metadata": {}, "outputs": [], "source": [ @@ -171,6 +181,7 @@ }, { "cell_type": "markdown", + "id": "45f0307a", "metadata": {}, "source": [ "Make an anesthetic plot (could also use getdist)" @@ -178,22 +189,10 @@ }, { "cell_type": "code", - "execution_count": 8, - "metadata": {}, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAELCAYAAAD6AKALAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOydeXxU1d3/32dmMklmsq8QyEYSgkERIbIjIrigVLoIfWx/hdZa1D5Vqv662PJra7W1m22pXZQuFh+fxxZsrTwoakFUFlkCBpQIJEMCSUjCZGayzZLZ7u+PO/dmJjthkoztfF4vXkPmnjn3nHPPPd/9+xWSJBFFFFFEEUUU4YJmvAcQRRRRRBHFvxaihCWKKKKIIoqwIkpYoogiiiiiCCuihCWKKKKIIoqwIkpYoogiiiiiCCt04z2A8UBGRoZUUFAw3sOICBw7doxZs2aN9zAiAnV1dUT3hYzoWvQguhY9OHr0aKskSZlDtfu3JCwFBQVUVFSM9zAiAkajMboWAZSXl0fXIoDoWvQguhY9EEKcG067iFeFCSFuEUKcFkLUCCG+2c/1VUKIE0KISiFEhRBi0XiMM4oooogiChkRLbEIIbTAb4AbgQbgiBBiuyRJVUHNdgPbJUmShBAzgK3AtLEf7fDg80s02pyYzF00d7gwxupYWppJYlzMeA8tiiiiiCIsiGjCAswBaiRJOgsghPgLsApQCYskSV1B7Y1AxKUS6Pb6+OO+Wl59v4kzLV24vf6Q6+lGPU9/bjbXFqSN0wgjH1a7m20V9awuzyXNqB/v4YwL/l3W4N9lnpeLSF6nSFeFTQLqg/5uCHwXAiHEJ4QQp4BXgLvGaGzDgs8vcd/zx/jJa6cxxOj4/IICfvypq3jx3vns/+YNbLt3PsnxMdz17BHOWxzjPdxRh9Xu5pm3TVjt7kv63baKep7YeYptFfVDN/6IYbhr8q+8BgqsdjcPb638l5/nSBG8VyJ5P0S6xCL6+a6PRCJJ0kvAS0KI64DHgOV9OhJiPbAeIC8vL8zDHBj/9W4db566yKO3T2fdgoI+1yelxLPlrjnc+qu9fPsf7/NfX5w7ZmMbDygvA8A9S4qG/bvV5bkhn/9KGO6a/CuvgYJtFfXsOW1maWnmv/Q8R4rgvRLJ+yHSJZYGIHjVJgMXBmosSdI7QJEQIqOfa5slSSqXJKk8M3NIb7mwwOXx8es9NSwsTmft/PwB2+WmGdiwrIS91a3sq24dk7GNF1aX5/LIimmsLs8dlFM3mbv4wrOHMZllTWeaUc89S4oiTuS/FAw03+A1GahtJKs9wgWr3Y3D7WPDsmKeXDOz33kq62Iyd41I8o1UBD9vq93NL/55hl/88/SgeyWS34lIl1iOACVCiEKgEfgP4DPBDYQQxYApYLyfBegBy5iPtB+89F4jrV1ufrW0GCH6E7568Ln5+Wx+5yy/33uWRSV96OK/DJSXAeCZt00DcuqP76hiz2kzUMWzX5gz1sMcFQwkmQSvyUBtRyrpfZSwraKeTbureWTFtAEPS2UdDp61BPbHv8Z6BD9fgE27qwEw6HVD7pVIREQTFkmSvEKIrwCvA1rgT5IknRRC3Bu4/jTwKWCtEMIDOIFPSxFSC2BrRT3TJiQyf0r6kG1jdVrunJPHpt3VnLPYyU83jsEIxwcK9728LBuH24fD7cVqd4ccJhtXluHxnaQkK6HPtUu5RyRw+CZzF4/vqOL+G0rYsKyk3/kGI5hz763uiES1x6XCZO7iO//4gOmTkrk3iONW5ra8LJtn3jb1++yC28yb0vKRW4+B9mXwvLZXXmD94ikAQ+6VS7nHWCLSVWFIkvSqJElTJUkqkiTpB4Hvng4QFSRJ+rEkSdMlSZopSdJ8SZL2je+IZTTYHLx3vo2PXZ0zpLSi4DNz89AIePFowyiPbvTQn7qn93cKd7arqgWDXsum3TWqAVJpCxCjFWzeW8u2ivph9RuMSDJsKtLXU29Wh8y3P3XXM2+b+OWuM2zaXY3T41cPhv7UHr3VhR8VPL6jiv0mC5vfOcuWA3Wqaks5DHdVtfT77IIPzFSDHofbx5YDtX2e/0gdREYL/RnclXkrqq8tB+pwuL389fB5Nu2uJl6vIT1Bz6bdNSFth4NI2PsRLbF8lLHz/WYAVs6YOOzfZCfFMW9KOq+838RDN04dNkGKJPSnslG+21pRz+a15awuz1Ulldtnyk5+CsemtP3lrjM4PX6KMo0sL8setN/g7xSMNYcfLIXtqmpRP1eX57JxZRlQxcaVZaQaerjzLQdq2bS7Bofby4M3lqrzyU8zAPD6B80hHH1vBKsLn1wzc1hcaiRwsxtXluH2yhILSH1UW72fncncxXdf/gCPz8+hWhuWLjdvVDVTF/Ci7K0uClaXDWSrGStY7W4eeOEY+2osONxe1i0oBGRJ5Imdp2i0Odnx/gWsdg8Ai4oVNbhgeVk2B89acLq9qmpsOGqw3us3Hs88SlhGCTveb+KqScmXrNK6bcZEvv3SB3zY1ElZTtIojW700N+mdrh9FKQbMJntPL5DtpkY9Fqe2Hmqz6GwvCxbJSo6jcBktqsHdHC//d0rGGOti+6t++9tAwi2Ey0vy+bhrZXkBgiI4vyozOPagjS+9uJxTGY72yrq+52H1e6mJDsRt9fPxpVlw7bBRIKtpigzgf/+0jwAlQt3evzMmJwcYpRW8PiOKvbVyGbTgnQDxxtsKlFZVJze5/mvLs9V13+g9RtNBB/k2yrq1bErEug9S4owmbs40dDOntMXsdo9pBpiWDs/n9tnTlL3u+IhN2NySh/njsHQe/3G45lHCcsooN7q4Hh9G99ccekJAG6ZPoHvvHySV96/8JEkLL3x9Fs1bN5bS9mEBPLTDNx/Q4lKbNYvLsR0sYvrf7qHn6+Zyaz8VHZVtahExeuXyE2N5/WTTVjs7j7ceyQZMoOJAsD9N5Qwb0q66v329Fs1VDV18uiq6aqkUZBuYMOyYtYtKAiReLZXNrKwKIOMhFhMFzv55G/3E6PV8MNPXkVRZgIgHxab3znLhmUlqoQUPI6hxhkJtglZBVTLwbMWDtXaWL+4MGQNnG4/bU43HzZ1kBirpbPbR53FwXVTM5EkEELw6KorSTPq+3DlipSorMtYzum+549yqNaK6WIXNeYuDHoNDref3R8288bJZn6+ZiZH6qzsOW1m7bx89pta+ekdVzMrPxWA1HK9ug7AZUsa4/HMo4RlFPDq+00A3HbV8NVgCtITYplTkMbrJ1v42s1jl5kmXOJyb+6oqqkTgKpm2Q7w+CtVWO1u6iwOlpZmqlz9hr+8x8tfWURjm5OkOB0dLi+5qfHU25zU25wcO99O1YV2ZuensW5BwbDHOFZqAIXIPfO2SZ3TxpVlbDlQy9FzNpVr/e7LJ3noxqkcb2ijzuLA6fbz8NZKSrIS2Ly3lq0V9ZjMdrXfQ7U99/jin4/w9y8vBGRVyvrFUzh6zhpQs/gw6LXYHIPPd7yJseLIsHFlGbuqWti0u0a9VtXUyb6a2hBpLxgKs5Fq0DO/KJ1Nu2v4zj8+4KnPzFL3ncPtxaDX4XB72XPaTEl2PdUtnWxcWaYSZRidfXHsnI27nzuiqrX+UdmI29fjR2Qyy1LWQ1srWTVzEhuWFbNkahb1NgfJhp6UTj1zkZ/p5Y73cp75SO8bJSyjgFfeb2LG5OQgVcel4ZYrJ/Dd7SepudhFcVbC0D8IA0YiLve36XpzR4+ums4X/3xEVV0cO98GQKxWUHWhgwxjDK12D5mJMpf23Ls9yVNdbh/ZibG0dHZj1GvZV2NhX42Fvx9r4M93zaEoM2HIjT9aaoDeNhXl/sFqGKhSD0i9VuD2SXh8Pt4+Y8Zq91CUaeR4g41DtTYq69uYlZfCsfNtJMZp6XT5MOo12N1+DDECh0eizuLggRfeo2xiIpv31lKQblAJdH+2ikiQ5oIJSVFmgiqt1bYe4aaybNn7SQCShMsjpzq6/4YSjp6z0eHyhvTl9UukxMewbkEBv/znGQD2myw88MIxHl11JQAOt48ndp5ibmEqcwvT2Pl+E/U2J73d1sO1L4Ln97UXj2O1e9AISIrT0ub09Wkfp9NgiBFs2l3N+sVTVJWnMr5gj0Dlme6tbkWSJPabLJc93kvFSNcpSljCjPMWByca2vnWrSOXNm6ans13t5/k9ZPNFGcVh3F0A2Mk4nJ/m643d1SUmcDfv7yQp982UVFnpcHq4GKXm26fREtnN4YYDcnxOnJTDTTaHJRNSORUSyd+CcxBXjB2d89LWm9z8t2XT/L83XPZcqCOTburcbh9PHjjVLVN8MF/qfMaDhTD+ztnzCEvfJpRrxrSry1Io+KclU6XD7dPPhQP1doAQX6abHNyun3E6TTYHB5cTR0AxOu0dOJTU0zotFrwyIfsvppWOlwyR1xncZAUp6Pd6WbJ1CwMep2qigtWAY2l8bY3wd1b3cq+mlYUB4OSrARqLnZRZ3GweW8tG5YV43T72HGiiQvtLgDKJiZRnGXk2Pl2td+iTCMms524GA0/2vkhOz9oUq/tq7Hwid/uoyA9gZmTU1hYlK4+E4A0Ywz331ASMk7ZgcSLw+0bkTu7MleFabLYK9U96pfol6jE6gQur5+qZlkiPd5gw2S2o9MIkuK0qnfYpt3V5KcZ+MWnZwbsMDKjMB7ZCEaqRosSljBjx/tyYoBbR6AGUzAxOZ6rc1N442Qz/7l0bAjLSMTl4E030OGlxCw4PT5VWgmGw+MHj5+Xjzf1udYfBHJOn5zkWL7zjw/4R2Vj4Epo6NLoGyxlg/v0SclcNzX0hQ9Wi3W6eg6YJVMzOHa+jUO11kC7GPUw1QjZuKvXCi52yQTV4ZY5+A6Xl8Q4LX6/hN3t53Rzh9pnh8vLsfPt/Gjnh8wvyuDJN06z32Rh3pQWipb02GNGU2oLfubKwagQlPWLC4nRCtXBQJG0QDa8g+xSHozjDTYeuGEq9zxfQbdXfq6xWsHklHga2pxsrQh1x4+P0dDh9HGioZ0TDe2sv24K03OSON7QDkgcqrVxpM6q2jBAfkYGvS7gQKK95HUxmbtY+8dDNLbJz6+2tSvkWfeGBtS5KCjNTuJMSxc2h4eXjzdTlpOCso/PWR089WY1T66ZyZYDdYDEugWFA2YjGC3GYaRqtChhCTNeOdHENXkpTE4dmRpMwc3Ts/nJa6e50OYkJyU+TKMLD3pLA9A3h5Fyff1zFSE2g6lZBkytTnx++QXSAKG5niHVEIPNIXPlGcYY2p0eAloSlXwcruvxDJLVI4UhfYy2wXLdggIMeu2AL7PV7qbR5mRSShw+v0RzRzfHzrcF1DJgjNEyOdXApGQ/521O2p2yRBKskwdUtVjwodXtldQ1UuxRTe0u1SV1UXF6CCc+WmvRP8GSx182MZHFJRmh61OGakMqyjTy6KorSTXoOXi2NSDJySidkMRjr1SFHMRVzV3Eafu632sE/ORTM3h8x4e0OWVJOD5Gg0Gv41CtlUXFGSGBpsGHsOLOe6kGfqvdzV3PHlGJCkBWYiydroGTyMbqBE6vxISkWHJS4jl2vo1Uo56PzcjhuYPnmJgUp47j4FkrkiSxcWUZaUZ9iCTeew5pRn1EePr1RlgIixDiM8DtgA+ZlftfSZJeCEffHyXUtto5eaGDjbddcdl93TJ9Aj957TRvnGzm8wsLh/7BGELZyHurzX3884PjM/ZWmzGZ7RSkG7iuJBOXx8vWo40hfRVmGqm3OkIO1A6nR/1/q90T0t6o11KUlcCJhnYmJMfR3O5iSmZfl+7RNlIP1L/y0jvcXp47KNuLZuUl0+31q0QFwO6ROWxjjBa7x6dKYsHQawVaTehhmp0Yy6TUOIqzEmm0OXn4plJ+/s8z7KtpZW5hGldPTqaqqYNNu6tVTny01qI/grVuQSEGvU5VhQWviaXLjclsJyc5DpPZzuf/dJjHAraRqVlGGtudGGN0VJ6XVUSxOhFCXFy+vgk1/BL8cnc1BZlGDtV2B2KA5BiQd86Y2VfTStnEJJXrB0IcBvacNodId0NB9vqq4Jw1lIgohvmB4AzMo7mjm6zEWArSDSyZmklBhpGzrXb21bSyvfICBr2WQ7VWHlkxLcTZIBi9CUkkefopCJfEskSSpP9Q/hBC/AYIC2ERQtwCbEJO6fIHSZJ+1Ou6CFy/FXAAn5ck6Vg47n2p2HH88tVgCqZkJlCSlcBrEUhYlA1s6XKzr8aC0+NXOSiAo+dklVdtq525han88JMzSDXoWf7zt/r0FSzNKOjn/FBhd/s4b5V/c/uMiVRf7BqTeIWh1A2KEbckO1F1A16/eArHG2xcaHOpElif+XhkSaS/Kbt9Em5HqAG7pbMbi71btT9cV2dldn4K+2pamTclDYNex74ay5jo4/sjWMFqwODcV0/sPKWqwDq75bWotzl5aFul6kUFYO92q6pAvwQp8f0bwYNhMtu52OEiJzmOZEMMm3ZXU1FnZXpOMvtNFqqa2lWvvA3LSvrEhFzKOn3r78dDpKuBoHiwBSNOJ3B5JU40yqrMr714nG33LlCfH0iDprUZyG443p5+/SFchCVWCHEbcu2UyUBYdDfDrCC5AigJ/JsL/C7wOaaQJIkXjzUwtzAtbKqrm6dP4Ldv1YzYuDhaUDay1e4mPUGvRhEr2FfTSnyMhsY2F41tLrZXNlJRZws5QAaDRsiHSn8QQJvDS0G6gQVFGVQ1dbJ2Xr6aUwkYFX3zQOoG5WV/89RFDtVasXd7eWTFNJVjnzclg027q1XVlVYjVDWgIqX0J630hzgtuHzg9YNWwGfn5mPp6gYh1NxSyqGjZCsYrwj74APyr4frVQ82xQgPsprT5fVTlGmg0ebE5ZWI02lwBQrheXzSkERFQWe3j85un2qz2m+yMH1SMktLM7n/hhJm55uR7RShruqXciBb7W5eP3lxyHaTUuLweP0qgQRIMej45ZprePPURd6oaqbL7cFktnPf80cBWL+4ULWhBBPm4OwBkajyGgjhIixfBj4JXIVMXL4Spn6HrCAZ+Pu5QOLJg0KIFCHEREmShmcNDhMO1Vo5Z3GwYVnJ0I2HiVuunMCv99Sw68MW1kSQmAuhHLwCh9vH7TNz2P1hC4frbExIiiUvzYDT4w/x0hkMApmoKAeMVoRKMMp/6ywOHtxaic3h4ZzFrqqZFIMshPflC1Y3KClGyiYmE6/XsGl3DXMLZcOwEIJrC9JUb6H1i6eokgsI1XAfPJdgotKfzUm9ptGg9fvxSfKanG7pVPtTYoKUTAaDZY4eCwQfkJv3nlXHuHFlGdsrGwHBnw/U0tXtw252qGvg8g40+8GhMCN6rSA7KY75RelUnrdxuM5GSVYi6Ql6lciO1Ni95UDtsBgAjZAdMGK1Grp9foyxGtocXh57pYp0o57mjm5Adt5Qnl+DzUG8XqsSl+Vl2WwNRN4/vLWSJ9fMjEiV10AYEWERQhQA/wkUAVagEtmucm6Qn40E/VWQ7C2NDFRlMoSwjHahr60V9STG6lhx5eWrwRRMz0liUko8r3/QHFGERanyFxwvEexho9fJQV0xWg2H62zMzE1hw7Ji/ri3lq4gt+HkeJ1qtFagvLh5afGcszox6jVYA+qgYG5WCZ4E2dhfb3PS1ObinNXOrLwULF3usEl6vQ+ih7dWqjE16xdPUbnip96sZs9pM1978bjqWHDyQjtev5zjSiNkDyanp+/hadALHG6pD1HJSY5TuXCXx4+fHgmnwSbfY1ZeMiVZiWpKFIiMCHuTuYu91WbWzs8n1dDjYGHQ61hdnktdaxcvH2/q97AOXqehJLr4GA2fn5/Pn989R2KsjnqbE/05myoZvX6yWbWJXF4JAjGs8STG6QJqL3n89m4/SXE61b18xqQkTK1dWO0ejHotdrePxjYXm3bXcKKhnSfXzGRXVYvq5BBMXCJdUlEw0uzGLwOn6FFTXQ28I4T4tRAiNlyDY3gVJIdbZXLUCn11uDy8+n4TH5uZQ3wgUjYcEEJw8/QJ7K1ppavbO/QPxgi9q/xZ7W4sXd0sKs5geVk2j66aztLSTJaWZgGw84Nm6lrtfR5Kb6ICkGHUkxSn48xFO91eP1aHF13AgJ0SFJ3s80usnZfHhmUlzCuSE/f988MWDtXaOHa+jc17z/Lw1sqwZLjtnS1248oyFhWns37xFOL1WvacNnOkzsqTa2byyIppfHVZCYlxWhJjtew3WbgQ8B7yS/QhKka9hlidBodbIilOR4YxJuRaUZBjgvJLCZkoK15J5k43m/eexaDXhWRDlh0p6votGKVgNDMBKzm+3jljxun2seVAHU8HJKlf/vMMh+oGtlUEr1OIRBf0tifGaTHGaHF6/PzPkXqcHln9lGqIoSQrgdzUeGblpXDO6mBhUY+n3PKybJaWZl6yN9i6BQUsKk4flKhokDMIuLzK85Sfhz3w/l5od9Hh8mLvlueXHN/zvCelxLHntJn7/+cYy8uy2bCshIXFGRSkG1Q7IkRe9ub+MFJVmFaSpD8CCCGskiR9SQihAx4ENgPrwjS+4VSQvKQqk6OBHcebcHn8oyJV3Dw9mz/tr2XPqYt87OqcsPc/HPTm2IO54TSjPqDukGMRdlW1cM+SIjauLONbf3+fpDidmpZlOGjt52VRjKCKCgHkF/Tv7zXy6ztn8bu3TWqMQ36a7G1ztjV8Rv3e3H9RZgK/unOWakgNdjteXZ7L6qcPqO7BsVpBvc2JDuhNRmO0Aru75wDtHWlud/vZWxMU6GeIodvnw97tp93pRaeR7S31NmdIMsZgzzTFA6p3sk8Fo6m3v/+GEjV1jbI/clNl++PLxxtVxqK3unMwBNveOl0+YgMuyDEawYSkWDpdXmwOD6+dlD3SclLiWL+4kDeqWtgf8JSDS/cGA5lYz85PUx0B+h1f0P87XF5iAnxm8Pwcbi9lExJoc3pZVJzBOaudQ7U2JqcaaGxzsd9k4fEdVcyYnKxmoghm4nprCyIRIyUsu4QQX5Ek6dcEGApJkrzAT4UQZ8I2umFUkAS2A18J2F/mAu1jbV/565HzlGYncvXk5LD3XV6QxoSkOLYdbRg3wjLQ4aPkpbq2II2FRelMn9SjivnuyydD7AlpBp2q0tIJCI4VW1yczqFaC0a9Dls/UkxvKKqIrm5fiFfR0tJM1dDZnw1opBhOhUfoURGazHbVVtIdOFGCZ6WM36jXDGqc7q1y8fgkldMVyERF8TSanZ+qSivK2NZfN4VFxemUTUwecB1GU2V2pM6K1e5hUXEGbq+Pw3U26m0y8VdUU3oNuP2QadSHZFoYCAoxVdDtk0gLpAUKRlaCnotdbvLTjLxR1UKdxUFBumHE3mAK1i0o4K9HzocwOYPB08/jvdjpJj5Gx4V2F1uPNrCoOIP1i6fQ5nBzoc1JRoJezWos22x7giOVXHTjEYV/KRgpYXkIeEQIUQHkBOwXDmA+YSwLPMwKkq8iuxrXBMbwhXDdfzh477yN4w3tfO9jZaNSP0WrEdw5J49f7DpDbaudwoyxryzZ+/DpnSJ+aWkm+00W9DqNSmymZBjYVwPX5CZz3dSsQAZba782hsO1Vtw+8HeHvoUKJ5uVoCcjIUZNhZEcr6PN6SUhVsvXbirlmXfOMqcglYkpPUGpo+2C2d+BrKgI04wxWO0eUuJjaHN6VHvQNbnJ1Jh7IrQ7goIe+9PbS4BWAz6//NkZpA6VkANDH1s1nb9WNODsJyDS4fayr8bC7Py0AY3Vo7lOvddoy4FaQOB0+1SDfmZSHI1tLrT9BD/2h962/Wtyk9HrtExIiuXNUxfpDOyh2ICocM7qUO1dN5VNUOc/0jmnGfXcPnMSm985O6z28TEQFJalMgIzc5P5+DWTOHjWoroaK5KQViPUrNcDVdAcz3o6w8GIbCySJPkD1RyvQzaITwBmAx8gu/+GDcOoIClJkvSfgetXSZJUEc77D4UtB+pIiNXxqdmTR+0ed87JRacRPH8w3L4Rw0Pv6oWry3N5ZMU0Nq4sUz8Vr6T1z1XwxM5TxMXoWFqayf9bOR2DXss3bplGUaYRp8dPr5g/9e+EWC2xup4tqRSYjtFpiI3p4YGyk+JYWppJV7ePV96XCz5daO9m0+7qMaua119Fx9XluSwqzlC59L99eQEblpWwtDSTRcUZ/GzNTEqyEtX2fqln7v1pgjTIRIWgzwlJseQkxwGyc0eNWQ6uUyptAtS12tlaUU9xZgJLSzNxBtzBx7qiYPAayRHkpTx441Q+PSdXVYk1trnU+JbBoKxTsA2qbGIiFzu7OVRr5a0zZpWoJMZqqbc5WVqaqWoRFhalc+/1RWGxT9y7pIjE2OHZUp0euYbMxCT5mcUE9ne9zYlBr1XHl5Mcz6SUOAx6Ech67Rs0O3UkExW4zNLEkiQ5JEnaLknS9yVJelCSpN9JktQ3IdS/KC52unjl/SbumD2ZxLiYoX8wQmQlxXHzlRPYVlGP0z2w6mSsoGzuoswE9fPJNTNZWpqJyWwnzRiDy+Njz2kzP9r5IU/sPMXjr1SxbFoWBemGPjEqSnR1m9NLt7eH8CjMaWObi/fq21V9eqfLy5NrZrJ+8RTqWmUppmxi4iUVQxoNyDr4FECOa0o16DHotTx38Dz7alr56+HzmDtlFcqsvBTmFqb1G6+j0F0/stQGYIyVX9WclHgutLsoyjQGXLgl1s7LpyDdQGl2Is+8beLBv8rquO/+70lZFy9EiLE63Mbf4fantPvrkXrV5pabGk+dxUFzRzc5yXFMSOrx/VEOJ51GqOukeBwCnLc6VAcG2eYU2B/dPtW1OV6vY8OyYp76zKyQWJDLIbJpRj1b7pobwgQNhjqLg4kpcSwsSmfFdNlrNEYrZFWqECwqzuDtM2Ya21w43PJElXITA2GgstaRYtCP+Jr3kYzn3z2Hxyexdn7+qN9r7bx8OlxeXnqvcejG4wAlq29RphGr3cPplk6WlmaqhnfZU6uWpDgtyfGy9BGnkw8CPz0caW5qPHcESX/ZibHMykth/eJCVgXKGC+/Qj4gd59qoaHNSVGmkXuvL44ITm7dgkJVNaionzYsK2bDshKqmjqpt8nj/enqq1HkFMWRMCHABSu0RlEFxuk0qm3lYpqmEMsAACAASURBVEc3G5aVsHltOY+smMa6BYXU22R1z2OvVPHEzlMsKc2Uc3F9bDqLijOoPG9jz2mzmmIl3DXRB+qv92GntDte38N7rrhyYiAZpZxtIth24Ud2JVb2UHZCLMHaZqNex5yCnsSScq2WGNbOy1NdduXUNro+0vblMiCz8lN5dcNilpZmsurqCUO2P3a+jfKCNCamxLNhWTE//OQMNiwr4WRjO/tqWmnp7Jl3QbqBR1dNH7Q/ZS0Vz8dIqHMfjGgSyhGizeHm2f113Dw9mykD5PQJJ+YUpnHVpGR+v/csn742t08OqUjBsmnZTExup2xiEpv31rJ+cSHxMVo1u/GJRpkTW1ScQU5yHDs/aKKz28e07ASqmrtYWppFqlGv1iZp6eymLCeJe6+XszwXZSWoSS4VP//Na8vHnaAoCE6br+jBH7yxFJALQTW1O/npHVezvfKCmhrE7QsNHmxqd3G41srDN07lyX+eUW0EcToNDW0yp1+UmaB6NCnVEu+/oYS3z1wEBNvuXRAoi9sKhKZcD7fBfqD+BsppZemS1VeKesrmcPP4Dnn8DQEiqSTfdHr8qrHf7ffT0uZhQlIsTo+PJVMzuOf6YrZXXsDp9hKv16l2CbmuiZcNy0r6zTwdDhRlJqg1VLKTDUG1X2Qouc4S47TceW0eILFpd41a9RPkGiu5qfG4vX5aOruZU5DK058bej/3Lr8cCXFLwYgSlhHi93vP0tnt7ZN5dLQghOC+64v48n8f47UPmrltRvgCMcMFOS36WZUjTE+IVQ9XJX1+UVYCqQY9yktWkG7gzrkT+PS1ueyqasHh9rFpdzVzC+W6IsnxOvXlUbLRXluQphZDUmqEp5brI4q4KClvgvM9HamzYjLbOVJnRZFLJqXEMSklno0rywKqM52qXlxYkkmNuUtdpzkFqWw92si7ptYQw66iitxyoE6tVqm4QDvcPnqnXA+3wT64v2BvvIEOu0/PyQvEe8nM0a6qloD7bzo/XzOTh7ZWMqcglRRjbMDQJqhq6mBfTSsF6Qayk2I5VGtj69FGirIS+9TheeZtU2Af1fBIoDx477xb4USaUc+3br2Ce5cU8cALx9hXY2FhUToP31SqFvKK1+twenwsLErH6fayaXc1s/JSSDHoVGI0tzCNeVPSQvru7d0YzLD0ZmAiyfU4SlhGgLpWO7/fW8vHrs5h2oSxq0t/8/QJFGYY+d3bNdx61YRR8UK7HPT2WAk+bB7fUcV+k4Xrpmaqh+6OE02YzHaqWzpJNcg5x2wBw/eUDAOHaq20O70qt63475+3OjCZ7TyyYhq7qloiNn/SQBy7UtN9UXEG+2pa1XxqSvaC4CqQt8+cxImGdu6/oYSfB6omHq6z9YnP2VZRr6bNV9arv5Tro43ecw4eY3BRtuBMDcH7ZsuBWuosstSyYVmx2m5RcXqgqqKcnXhRcTqz81MHlJI2LCtWGZyxyrGVZtSr8U3KuFbOyMHmcPNyZaMqeZYXpIWU5VbyyOk0IhB3JNR1CR47EDKPSCMmwYgSlkuE3y/x7X+8j16rCUt6/EuBViO4d8kUvvG393n9ZDO3hDF9TDjQ30YPDugKVsekGfVsXluulnWVD8aedOaz8+XswFVN7WpdimCVz5E662XHJIw2enPswfmzNu2uYf3iQlX1AyKE8Myb0qJGzu85bcbjk9hX08qi4gxm56f0me9A0slYY3CVjCylOd1eQFLVVKH7JphZEiEqn8Ul8v5RUsIEz7G/zL8DlcseTQTPRX7O1RRlGtU4mlUzc7h9phyPlptmwHSxi4dvKuVInTXAcFygos6q5tbrb+yRuNd7I0pYLhFPv2Nif42FH37iKrIDLoRjiU/Nmswf99Xyw1dPcX1pFnExw3N7HC8Ep38JDl5U6mMo36WWyxKL0+0PJOMrCNgILOyqkiOkFZ02EFINMFK5tqFURA63T61Zr6i2lPY9EeGDFM7qhcEKj40VBuOilVotirpzw7LiPvE16xYUoOR8VtZkIJVP8JoOJpWMF2evPOdrC9J46s1qNq4soygzQSU4j6yYxvcD9Whm5aditbs50dDGflP/ZQ8iWULpjShhuQRU1Fl58o0z3DZjInfOGR+uQafV8L2PTeczfzjED175kMc+fuW4jGO46C+gSykEBj2pRoKN3P399qOO/lREVrt7yCqUINcQ6S9YbrD+IxHBtieDXovD7esz5v72Qe8Dtb+UNZG4V4LHrTBE0HeswfPpzYSNd5bqkSJi3Y2FEGlCiH8KIaoDn6n9tIkTQhwWQhwXQpwUQjw6WuNp7ermgRfeY1JKPE988qpxtW8sKM7g7kWF/NfBc7xcGZnuxwr6D+iS1y44v9WlINJ89ocDxc1VKeCkRMkP5iKtqAcNeq0q6Q0073C50Y4FlHmvW1Aw4JgHm6tCRJ0evxqf81EJHIS+70QPUyB4ZMW0EMneYnezsCj9khNmjjciWWL5JrBbkqQfCSG+Gfj7G73adAM3SJLUJYSIAfYJIXZKknQwnANxe/3c9/xRLHY3L967gKRRDIYcLr5+yzSON7Txf7cdJ8WgZ8nU8GZsHk0MVC9+uLW8B8rTNZ6FrYaCkmrlUhIIDpRKp/dvI33uA2Ew1c5gElhwyhp5LavUw7g/RPr6DJSmZVtFvZo6ZnvlhVF3xAjnOkUyYVkFXB/4/xbgLXoRlkBxr67AnzGBf8PMkzo8eH1+vvG3Exyps/HUnddw1SgkmhwJ9DoNf1h3LXduPsiXnqvgh5+4ik/NmhRxnmL9YaADZbi1vAfK0xXpKoPe5QaGQu91Gm68yL8CBlNtBavUTjS0D5nFOtLXZ6D3YXV5LnurzYEcYmE91vpFWNdJkqSI/Ae09frbNkA7LXKhsS7gx4P0tx6oACry8vKk4cDc6ZLu3nJEyv/GDunXb1YP6zdjDZu9W1r99AEp/xs7pPXPHZFaO12X9HuDwTBKI7t0WLq6paffqpEsXd1j+lsFs2fPHvFvh4NwjHGs+h3ttQgXhjP3y12f8VyL0dozI70XUCEN4/wWkjT6lHAgCCF2ISew7I1vA1skSUoJamuTJKmPnSXoegrwEnC/JEkfDHbf8vJyqaJi4FyVlq5uXjh8nmfePovT42PjbVfw+YWFQ01n3ODzS/xh71mefOMMiXE6/s+8fP5jTi4Tk+OH/K3RaMRut4/BKCMf5eXlDLYv/p0QXYseRNeiB0KIo5IklQ/VblxVYZIkLR/omhCiRaldL4SYCFwcoq82IcRbwC3IWZZHjN/vreXpt00sm5bFt267gqIxSNlyOdBqBPcsKWJJaSY/fPUUv3qzmnfPWth6z/zxHloUUUTxb4hItrFsR65E+aPA58u9GwghMgFPgKjEA8uBH1/uje9aWMAdsydRHJTi/KOAaROSeO6uOdRbHbQ7PUP/IIoooohiFDCuqrDBIIRIB7YCecB5YLUkSVYhRA7wB0mSbhVCzEA27GuRXae3SpL0/aH61mg0Unz80GqifwfEx8dTUFAw3sOICBw7doxZs2aN9zAiAnV1ddF9EUB0LXpw9OhRSZKkIcNUIpawKBBC3AJsQiYef5Ak6Ue9rq8CHkPOsu0FvipJ0r7B+hzKxvLvhKj+uAdRe1MPovuiB9G16MFHwsYyFIQQWuA3wI1AA3BECLFdkqSqoGa7ge2SJEkBCWYrMG3sRxtFFFFEEQVEcOR9AHOAGkmSzkqS5Ab+ghzfokKSpC6pR+wyMhYO31GMGprbXXzxz0cof3wX39t+Eo/PP/SPoogiiohCpBOWSUBwSbSGwHchEEJ8QghxCngFuKu/joQQ64UQFUKICrPZPCqDjeLS4Pb6qbc6cHtl4nG41srKp/bx7lkL1+Sl8OcDdXx3+8lxHmUUUURxqYhoVRihObQV9JFIJEl6CXhJCHEdsr2ljxuzJEmbgc0g21jCPM4oLhGV9W3cvaWC1q5utBpBTkocF9pc5KUZ+J8vzWVqdiI/fPVDNr9zlpUzJrKgKGO8hxxFFFEME5EusTQAwTkdJgMXBmosSdI7QJEQImJOIUmSaG5u5vfbXuU7P/01f/zbTjo7O8d7WOOKTpeHe//rKPF6DY9//EruXTKFsolJ3L2okJe/spCp2bKb90M3TiUnOY4n3zhDpDuZjCU+ikk4RwP/TuvwUZtrpEssR4ASIUQh0Aj8B/CZ4AZCiGLAFDDezwL0gGXMR9oLZrOZR3/2a7o6O8jOmURufiGz5i3k/NkaHnn8ZzjsXSSnpvLVuz9Hfn7+eA93TPHcu+do7nDx0pcXcE3egMkUiIvR8uWlxWz8xwe8e9byLyG1BBek2lXVon4qif+Gkwgw0nNfDYXLTXbYO23+wbOWPkkoIz3xZO99MNTzv5RnHglzj2jCIkmSVwjxFeB1ZHfjP0mSdFIIcW/g+tPAp4C1QggP4AQ+LY0je/vi7kPs/Mc24uIN3L7ms6RnZoVcn5xXwILrZU2dzdLK7/7rRVqaGtHpdEycnEt6ZjYLphdyzTXXoNNF9OMZEbw+P88fPMfC4vRBiYqCO2ZP5mdvnOZ/Dp3/yBGW/g4PpRbNO2fM7DdZQsoQ37OkaFgHSCTWHrkUDJadOnitlLYDHbIblpWoJX77K9UcqcTXanfzwAvvsa+mlTdPtXCo1obD7eXBG0sHHPelPPNImHvEn1ySJL0KvNrru6eD/v9jwhBtfzmQJInHNm3GdPpD8gqn8IWvPERc3NABmKnpGdzxubvUPpoa6rGYW3in8gybn9+GEILu7m6umTOfDZ9f85HIXDwU9la30tTu4nu3Tx9W+7gYLZ+4ZhLPHzyHpaub9ITYUR5h+KC84MHEQzEbTp+UzHVTM0PKEINclvjgWcug9Tc+SpUE+8Ng2alD14phHbLBNeZ7t1Hq30SS5CJXRm0N/CXvh6PnbFjt7gEJyKU880hgPCKesEQa/H4/lZWVvFlRxYWG83S02/C4PSxceiO3fmLNiPsVQpCTm0dObh4AS266FZAJzqG9e7j7ga+RkJTELx//fx9pArPzgyYSY3VcXzr8+jGrZ+fy7P46Xj/Zwmfm5o3i6MKL4MNtxuRGHG4ft8/MCalF01tnvquqRa0xsnFlGdsrL+B0e4nX6/pUkVRKPDvdPuL1mnGtdT8S2Bx9yzUHE1qbw90vke2vnMCWA7Uo5YwBVfp5fEcVe06bQ9RlY6kq6n0vk7mLvdWtrCmfTM3FLjw+P7PyUthXY2HLgToevHHqZdcYGinj8e9SjyWi4Pf7+dqjP+Zi8wWuuqacSbkFzJwzn5TUtFG9rxCCedfdwLzrbuB4xUHW3vsA06ZfzbcfuHtU7zsa8Pr8/LOqhRuuyCJWpx32766YmEhemoHXTzZHPGEJfjkVpBr0GPQ6nth5SiUqvWu1v3PGzPScZFweH7mp8ew5bcbefYLDdTa1nxMNbSG2hC0H6ti0u1q97nT7SE+IHTPufKQHUW/p5OBZC/ffUMLBsxauLeh5nxQiO29KC6nlepVYbK9sxGb3cLqlk9LsBM622gM1S+Ty4eUFaWzaXa32X5RpDFGXjbaqyGTu4rsvnyQnJY7DtVbqLA6eecfEH9Zey5NvnGa/yUKKQUebwwtAQboBkJ9fsHSljHNvdSuz81MGZBzCRRDCuS5RwjJM/PDXf2Rq2VXcede9w2rv9/v5wZ93YDl3Cme7FUnygxAwiPlHFxtHUtZk7vvkDeQVFvWRTK4un8fV5fN49+3drLvvq3zv6xsoLIzcdP69cbjWis3hYcWV/VVKGBhCCG6ens2fD9TR4fJERAXPgaC8nA63Vy1C5XD7AIkNy0pYXpbNAy8cY1+NBYfbxzW5KcTHaNhvsrDfFOpz0tTuAmBuYSoxWi17Tpt5eGslT66ZCcDRc9aQ9lVNneyrqQUGPxjG+yBaXZ6Lw+3F6fHj9vrZc9rMeasDk9lObaudOouD1082o9MI1l83JcQ29cLh89RZHGpfh2rlNShIN1BncbDfZKG8IJUNy4pxevzMmJzM7TMnhdhthqNuHAkUgnLOYqfe5lS/12kEVruHNc+8y9QsIwBtDi9zC9PQaQQP31TKkTorDrdXXU95jXzMLUxjX00r+2paMeh1gxbI68+J4VIQThValLAMAw6Hg+qqD7hl1R1Dtu12ufi/jz2Jo81MVtFV5F1zHfHJaWg0Q3Ponm4nHS31PPWXndgaTcQaEvnRxodJSk4JaTd/yTKunFnOj371NFPLruLhL312xHMbS+z8oJm4GA3XjaCM8k3TJ/D7vbXsr25lxVUTR2F0l4dg4zOAw+1Tq0U63V42761lw7JidlW1qNy10+3le/97EqfHT1KcjoJ0Axc7uzHqtTS0uai3OZlTkMq8KRksmZpJU7uTPafNLP/5WyyflsW+GgtzC1O5OjeV+BhNnwN0IISLM72cg0ghuouKM9iwrJhrclN57JUqrsxJps7i4Nj5NrXdginp7K2WbRJ1Fge5qfG0drpweiUMOg3/Z0EBn742l+2VFwCJdQsK2VZRz6bdp3hkxTRSDaEHbbAkVLRk5CUxFFWkzeHmg8Z2zrR00tXtA+Q4Dj+QkxzHIyum8eDW43j9ElXNXcyYlMS8KRkh6stZ+alY7W4Mep0qyW7aXc2i4nQAFhalD7jOq8tzVelssEqaQyGctrsoYRkGnn3pdeYvWTZkux9seYUz+3Zw5Y3/QVL25Eu+T0xsPOl5U0nPk2tbOzusfPWb3yFv5iK+d0+o/SYxOZkvffUb7N39Op/54n3cvOoO1t0+9BjHC5IksevDFq4rycSgv/RtNzM3BaNey35TZBGWHtdXn6qWWl2ey9NvmZhbmEpuqoHjDe0AON1+nG43OclxXGh3seXdOrIS4wCI1Wo40djRp//zVgeH66p57t06bA6Pyv2+eKwRgBitlnuXFGFzuPnOPz5g+qTkENtFf9xruDjToQ4i+eDtsX0o9o37nj/KoVoruanx7KtpJUYrcHr8mMx2mtudIX24vH7u2nIEf0DQ12sIkQYcXj8732/i09fmsm5BAU+/ZeKBF47x0I2lbFhWgsPtVVWGCkevSAMOtxer3T1il+e7txxRCWBvKImIPF4/7zd2sKQkg90Bp4TTLV3ExWg5XGfj6Lm2EDWXLJ3JdrP1i6fg8vqQJPj+x68ccJxpRj1PrpnZrxPDeCFKWIaBE0cPs/aeBwZt853fvkBrbRXzP/MQQhOeuNP4pDTmfPoBTrz6HI8/G8fGL9zep83iZTez4PrlvP7yi3z2S3/n2d/8Ar0+8gy4Jy900NTu4sEbp47o9zFaDXMK0zhQM+4hSkAPQbF0udm89yzrFxfyyIppKre5ee9ZAA7VyjaShUXpHG9oU1U3AN1eST0kzXY32Qmx2Jxu3L4edWlzRzcANoeH5Hgd7U5ZL++XZH+ifTWtPP22id0ftmAy29lvslDd0tnHYB2MsfAqs9rdPLy1UvXwUuxD2yrq1TXw+SUK0g3sOW2muqULALu7b244f5D22O2HGK3A45MQyGk46m1OvvX39wFJXe8YbTXzpqQH3JKL+7glG/TagM2rf/XSUNhWUT8gUQmG2e5W94KCbq+fw3U28tMMqprL6faTnqBXY3MAFhVnqN5jj++oGjJWJ5I8BSM98n7c0dLSgs/nJTYubsA2P932Ni1njnP1ys+HjagoEEIw49a11BzYid/ff0JGrVbLrZ/8NJ+4cx3ffPxnYb1/uLDrwxaEgBumZQ3deAAsLM7gbKudpl5c7XhAUSdVNcnSSHzggEoz6llelq0aZOcWprF2fj7nrY4QoqIgIbbHjtbS1a0SlRitoCjDQGKsrEKN02mYlBLqwq6ct1sO1GEy28lPM7B2Xj4lWYksKk5XD9LxiNreVlHPntNmclPjVWcERVU4Y1IScToNF9pdqmqroU1+pjoNxAZeoXht/96PnsAaBV8+e7FLJSop8TFsXFnG6vJcHlkxjXULCnlyzUyV8IMsrW1YVozD7bvkdbHa3TTanMTphvbOjNdBmiEGpWlWop6sBJk4zC1MY26h7KxwvMHGEztP4fT4mVuYxpyCVKZkyHtIIb7bKupD+lb2YO/vIwERT1iEELcIIU4LIWqEEN/s57oQQvwqcP1EIPo+LJAkia8+8l0+t/7+Qdudfmc7M25bG67b9oEQgqySGZhOfzhou4KiEqytkZlgc/eHF7kmN4WMy4hDUQIk90eA1LK8LJulpZk8dGNp4PAqUK/tqmqhzuKgIN3A1bkpnG7uDFHfBKOrWz4kjfrQV9HjkzB3uekM6OxdXj9nL4amAkoIEJ1ur580Ywwv/edCJqXGs3nvWWbnpwapgmrH/ABaXZ7L0tJM6m1O6m1OlpZmUpqdyOqnD3CisQOX168ezH5JIjleVp54/ajMmdM3eJyzN+iyOUAc9BrIS4vnmbdquO/5o1i63NS12nl4ayXLy7JVjl/+FGzaXc2WA3XDnpfJ3MWqX+/juYPncHmHjsN2esHq8OCVQK+Fi51uYnTy/D640E5Lh+ygUW+V7WlIEodqrRyus1FRZ6Ug3cD3PjY9hCgqUAhnpKi/gjEqqjAhxDcCgYuX289w6rGsAEoC/+YCvwt8Xjb+9NLrFE8rw5gweIlijUaDTj+6gXsarQ6/zzdkO11M5Gk3m9tdvN/YztdvKb2sfqZNSCTNqOddk4U7Zl+6DSuc2F55gT2nzcyYnBKi3jOZu9j9YQuTUuKoszjY/M5ZYgbgvIOhCXgAagUo52mHyxvSxhX0+CelxNHY5iIrQU+7y8ucgjS2HKjl9pmTsHR1c/RcG1MyDDx38DzrF08Z8wNI0fsH21hWP30Aq92DRsjqrZunZ3Okro3GNvlwNcRocHj8KtEZzsENPWsBsqrsRGOHaq86VGvlLxXn6HD6cHs/4L+/NC/ol1Kvz6Hx3ZdPDsgkDAWvH9bMnqyqtxpsTjpcXgRwod3FhXYX0yYmkWqIwebwUNUsqwe3vFvHs1+Y06e/SFN/BSMsp5AQYmvwn8BMwhMNr9ZjCdxHqccSTFhWAc8F0rgcFEKkCCEmSpLUdLk3v+sTN/OF/3ydupozFBQPZhsY/YBF6/kzFE37wqBtPG43zgisgLj7VAsAy6+4PPdOjUYwf0o6B0ytSJI0zoGifQ8lk7lLPTyD4enFeRv1mj62hM5uH5ogoqIRkJcaT5217yGWmaDnqTtn8dSb1aoN47WTLbx2sgWDXsfJCx3sN1k4a+5R3471AaTo/4NjL/7fbWXc999HcXrkuR+pa6Oru4d4KrYUnWb4RCUrQa8SZQWTUuJIM8RgdXhobHPR4ZQp8vRJySHt1i0oVL2whouc5JEzkH4Jdp1qwWr3kGaM4dGPTefrfzuB0+NXie1Zcxc2h4eCdANzCtM4Z3FQkpUwYieD8UK4VGEdkiStCfxbDewKU7/DqccyrJotI4EQgj8+9TP+8ufNeD2eQVqOfmoySZKGzB32X5ufYvXayAuc3FXVQm5aPCVZI3ftVDC/KJ2mgG5+vCDr5AUblhWzbkFPHNHjO6qw2j3E6TSsunoiMyYlkZWgJyshlglJPQeSksouJug81GlCjdR+CRoDcSy9Ye5y89oHzcyYnMzaefnMykthxqRE5hamcW1BGt5ARxOSZcISrx9+MOpIEGzDMZm7+MKzh3n6bRNP7DzFw1srsdrdsofYu3U4PX7V/nSh3UWHy0uqIYayCQm4AnV5vP6B36f+WIl6m1NVC87KS0Gn0fD+hU4mp8bLNq55eWxYVsK9vYirwvEP98C22t28daZ16IYDIEYr+PnqmSwtzcRq91BjtnP91EwE8vNeWprJo6uu5JEV0/j7lxfykzuu5oZpWWzeW9tHjRnp2Y7DpTf5Qa+/vx2mfodTj2VYNVuEEOuB9QB5ecOP3tZqtcyYdS0XGs6TVzg+YqfMnQ/OAxw58A4TJ+fx8evCZmIKCxxuL/tNFj47Ny8sEsaCItmv/4CplcIM42X3NxIoMQaPrJgWcihtXFnGeWsFJrOdDpdXVcekGmKwdfV4dTk88vbUaICAesvbj1+GNMgBW9XUzr4aCxuWFZMcH8Oe07KH0tdePI7JbGdpaSYbV5YNK67lchEcF6PEU7i9/hBPLEA15ifF6chKjOViZzez8lL4w7pruXvLkWHdK3hFjHqNSowMei3ZSTIhPWd1kGqIUY3586akj9gbMRhbDtRxsbN7xL836LXkphtU7ziH28vOk7I0Hx+jYePKsj4xNwO5hkdCosnBMCLCIoQoAP4TKAKsQKUQ4n8lSToHIElSX/eXkWE49ViGVbPlcgp92e1dGIwDH2J6QxKuzjbiElMGbHM56O5qJzYhecDrkiSxb/frPP/734zK/S8HbwcOmRsvUw2moDDDyISkOA6YLHx27viUGwh+2YNdPosyE9h27wK2VdRzbUEabq8fk7lLdRlud3pVV1kAzxBVl70SqookGIlxWnKS41m/eApHz7Wxr6aVWXkp1LbaVaKiuKYWLUnAZO7i4a2V6sEV7KIajij84PW4tiCN81YHD99USkGGkS0HanG4fSyZmklRphGTOTQq/cqcJGwON6eaZSKsBBb2hmJ/CYbbJ+FxyzaKDqeHi51upmYZSTXEkBawU8iQBkxTf2kYuWZCg/z8v/jnI9xUNoF4vRzQarN72H6ikTaHl2/9/X3qrQ4utLuwdLn51m1XhPQRPAeH28eGZcXDYhrGI43+SCWWl4FfAa8Bf0Je8a8JIXYAD0mSNHKyHooh67EA24GvBOwvc4H2cNhXgmFubiIja+A0JHnXLObs4X9Stmx1OG+rorHqMNklMwa8fubk+0y/elZEJqfccaKJdKOeOYXhyakmhGBBUTpvnzHj90toNGNg3+r1YgYbTZ8JqHwAVa2yujyXh7dWst9kITc11EU42N6SlxbPOYtz0OOqN1HRaqDT5WPr0QbiYzQ4PX4WFWcgSRI2h4eiTGOfeAclESNUOxFhhAAAIABJREFUqbEdynjDwfkGr8e2inpMZjtH6qwUZBg5es7GvhoLR89ZMZntzC1MBQS1rXYudnbz0nuNgbQ3MtEYiNb2JirQs5YSqDaZmot2/KASlZzkOJxuv5pGp3eZgkvBugWF/OO9C5yzDq2G7c0QJMbpQMiZA5S4FoNex/c/fiXrFhaw/rmKEHd0xY09OF3LjMnJav2ZPafNfSTmgTAe0s1ICYtWkqQ/AgghrJIkfUkIoQMeRJYK1oVjcMOsx/IqcCtQAziAwS3cI4AY4vD6+prruefAa9gaTKRODu+Dc3bYuFh9gp98/b4B29SZqvnUTYvCet9wwN7tZfepFu6YPRmdNnye7fOL0vn7e42cudjJtAlJYet3IPR+MYMjym+fmQOEqiq2HKhT05U8dONUvrjlSBD33IPspDjqLMP3MDLqtdjdPSlDFCP47PwUbp85icd3VHH/DSV9uNONK8tQsiUrqpbeKpZwqcuC+5PTw8uu4R6fnCtNicBf8/QBLnZ209nto7PbqQY7DoZYraC7HxdkvVao8T+9yc+E5Dj1IFfUg8FlCi4FaUY9f/rCtXz+T4eH9AzrzRC09/LwK0g3sLo8F5O5i/XPyerT3NR46m1OJqfE8+iqKwHZrX1rICZoxuQUNiwrwen2MWNy8rDnMB5p9EdKWHYJIb4iSdKvCewHSZK8wE+FEGfCNjqGVY9FQlbLjQq8Xi9ej3fIdr998gc88K3HuPBhBVcs/RSayyzS1WVppubATnweN7/5xU8GbWs68yET142OtHQ52H3qIi6Pn5UzcsLa73zFzlJjGRPC0l/9DyU62qDXcs+SohB1kzNw+Ht8fpINMayePZl/VF7gYmc3Oo1sT5mVl8LVk1NosDlVV1kFOo3o14BtjNVROiERc2c38XoNZ1rs5CTHqZ5Xz35hDr/452k27a5RC0dZ7W52VbWESDHBXGu4XVaD+1NSpxw9Z2VfjUUNjv3FP08zMTkerbCpXnDDUTKlGvWqWlFB2YREnvrsLP56+DwHz1qptXTR6fKph3RxphFLl5slpZl8dflUVT04UhRlJvDJWZPU5z8UdEDw6WGM1SKExD3XTeHhrZW0tMvJN7MTY1k6LQvTxS6+//ErKcqUPcEe31GlqjfXLSgY0L43GMbDLXmkp99DwCNCiAogJ2AYdwDziYCywOGETqdDqx3aq0ar1fKbH3+PH/33G1T8/XcYU7MpWXQb+vjhG5jdTjvnjr2NrcFEQno2T3z7oSHT8p84epj8KcUkJw9sgxkv/O/xC2QlxoakQg8HJqcayE83cMDUyl2LRj+7c3/1PxxuLyBUYqOom9zeD/D4ZL75UK1V/T4/TfaEmp0vJ41842Qzm/fWqqqydEMMloBU4xvAaO9we9U0IusXT2FSSicbV5b1OmBCC0eNp5E3zajnwRunhtgGgtO8KIjXyYGE/SE5Xkd8jJbmjm6yEmORJGgJMqB3uLxsr7xAvF7LvKJ0TjS2s6g4nUdXXcmuqhb2Vrdyzuqg3uoIm31h3YJCnB4/+6rNVDV19ttGiUdSppWTHIfPL6lj//Frp2lzeogNxDg5PT6ee/ccIMfKzM6X7bVKIlOFMYiEIl7DwYgIiyRJfuAHQohfAMuR41ZSgQ8In0dYxKBo6hX8zx9/x+rPfZGYIfJwffOzN8Fnb+Kn297m5Bt/wdvtJDFrEhmFZaRMzEenD00N42hrpenUMaznz6CLM1Awawk/+9bgkf4KXE4Hf/vvZ9nx4gsjnttooaXDxZunLnL3okK0o2AHuX5qJn85Uo+924sxdvSDQnvbWR68MTTYU1E3lWQnsvkdWfWyqDhDVb2UZify2CtV3LekmMdeqaLO4iDNGKNGpa+bX8CDWyuxOTwq964YrLUa8Pmhq1tOoz5vSprq5rytoh7KUA3S6xYUcKKhLSSFymikiL8UBNudFM+wepuTnOQ4Vl6dQ0WdlWPn29BqBD6/pKq8YjSywXtSchzNHd3otIJVMyepqq1UQww3XJGlJv9U3Jhn56dSlJlA0ZIElpdl892XT4Y1FiTNqOdbt14BXMGxczY2/OU9nG4vrXYPE5JiEQiaOkKl0AvtLpLi5H2aHK9jydQMXj7exA3TsviwuZOkOC0nGjspSO/JH7ZhWYka2Nqfo0Uk1LYfCJf1RkqS5EA2nm8Pz3AiE99+4G7+tucIv/3Z40ycnMdtn/z0kNH4X1u9BFYvQZIkmhsb2PTX12g4cQCfJ9jvXCI+KY3sqdfw46/fe0nG967ODn7708f5w69/EZFG+61H6vH5Je6cMzqFuVZcNZEt755jz+mLYVe19YehOP+izASe/cIcrHY38TEaFPuLcuArRu0t78p5vYoyjfz0jqs5UmdVD12bw0OKQUdcjI7mdhc6rYb18/K55cqJPPjXSs5ZHcRohar6UhwHehukFXfW4AqKl5si/nKh5A5bWppJSVYCm/fWsrp8Mg/eWMpbpy7y0LZKvrykiN++bVIDTBV7fUNAVXjsfDvl+WlsWFZMcLXIVIOeijor+02WgMqoR4otykxgcUkGT+w8RXpCbNiltln5qez9xg189vcHaTVZaO6QXaibOlyqTcwYq8He7afD5SXNGMMf1l7L22cuAjB1QiIz81J5YucpFhVnUPb/2TvvuKqr/oG/z71w4bL3FAHBhRtx5CjNUZZZWdrWllY2bNnee5dZPWo2tJ4s7WlYaZaWOweOXCmCgAgqyN53nd8fd3hBuFzgXrj24/16KeN77vec77mX8znnMyONql1jSv26FUPrfwZd2eXY9fJ/uChXjR7EVaMH8cPG3XyxcB7e3j5ced3N+PjZ1vELIYjsFMNrD81w2FgO/L2Lld99w0dvv0pISIjD7usodHoDX+/IYXhiMHFOijUZFBdEiI+KVftOtolgaY4KwrqmhrHol9HmMntMIpP6R1uMx+Y6HGA88WSe3kFWYRVTUyL46e88ymp0rD2Uj1rlRnZRFQmh3mxKL7QU+zKPxaxqNJ9KzKq7BeszLIt5e6tOzPYWkEzqH02wj4elHv3GI6cpqtRahEpcsBfjkyJIzTaeZMpqdHi6GWNW1Co3S0xKRkEFL/1sdEow2x8a2r23hfrohSt6c9vnxvcvJTbIFFtkFPaVtcbyw0WVGrIKq9iRVcSk/tHsPV7KBd3CWJ9WwMyRXTiQV8rCjUcZkRhiUYXZeg5XVot1CJZmcsXIAVwxcgC5ubk898ZcfP0DmHzDzXh6qpt+cSvJP5nHN59/TGyXRJZ+Oh+FgzMpO4qV+0+SW1LN0xOTnNaHUiG4uHcE3+48Tmm1Fn+1c6tK2msAtd5Fmv/gzanQZ49JZMWeXEBQXFXXs2zNwVOMTwpn4cZMcourqdYaE0tmFFQCkscn9KhzAjGnf7cWIPVPJdYLT3urSoK8VWelqjefuKYNjeV4cRX9OvlTWKnl+ct7WYzXd32ZyrbMYmp0ZyL2zdHmZm8qOMhntwxu9P1pC+N1QqgP380aXqcmynu/p/Hj37mUVusorNCQXVTFiMQQCiuMtXPMFUPN5ZPNJ9nGKkbWf45/fa6w/49ER0fz8dw3+H79Tua+/Ay33fswIWGO1WNLKTl65BA7Nm+gsCAfH18/PnjjRXx9bavh2hMpJR/9mU5imA/jnazXn5oSw5dbj7FiTy43nRfn1L6aon4FSesaGRkFFezMLmHr0UJLNLjZDmL9vVmnbg7iayiYr6GCTo3tXF1l4Wlobqy/Vml0ZBVWkVVYxeMTeliEyvLUHPrFBLIts9h4iukVwdy1R/AypagxL8RPmTYw7W1zqD/f0YFqSqt1FqExunuoJRYFzrg/w0GLe3pSpB9jehhfX6XRkVFQ0YqAzvajQ7C0kisvGMj4lLnMuO9hLp18Ld17Nx7I2BR6vZ68nGzSDx3kn317EEIQl9iNR2bdSlhYy+uYtCVr/snn0Mly3p7Sz+nBi32i/ekV5cdX23O4cWhsu9qabOm7V+zJs2S0HZ4QTEpcEJP6R9G3k/H0Mql/VJ3YCvNJJtDrzEJlvWjWv7+rCJDGaGxuzlRMzGTmyC6oVUpLNgNzQKN1ATWAYCvPKKh7GrPH5tCWwsc8TusNghFj1I7ZFtQ13BeNzkBSpC8LNx6tI4DMJZytn6m9Bag9dAgWB+Dt7c2XH3/Is299wO+//EBAUBCjLppIp85xDbavrChn/56dpB3cR2XFGXdFpdKNiKhOTBo9lAdvv94uN2dXQqs38Nqqf4gP8bYEDjoTIQQ3Do3l8e/2sSWjkOGJ7Wdvsq3vNvp5DU8IZt71xlxuy1NzLDXqA71UFpuMdQVBc4xM/WqMrixEGsLW3Jhjgqxdahesz7AEVqqt7FXWQtWsDrMuw2yPzaEtDd7WAt9aRWmdt2zB+gyLF2Hm6UpLUbS+nQIsp9f6AZ2ubLQ347KCRQgRBHwDxAFZwFQpZXG9Np7ABsAD47N8K6V8tm1HakShUPDiI8byxfn5+bz/2VJ+/PoLVCoVvv4B6PU6SouLEQqBt7cvvfol8/yc+wgIcE5+sfbgi7+yySioZNG0FNwdGGlvi8nJ0by3Jo0P/khvV8Fi69RgnZ69MW8uwFRGt6vF68k6INNVjPAtwdbcTEmJscyD2W5kbeg3G+XrL6TWqU6sBW5TC62rGbynpMSw8YjRppJnymZtDoY0n0bqe/O52jM0iJTSJf8BbwCPmb5/DHi9gTYC8DF97w5sA4Y2de+BAwfKtqKmpkYWFhbK4uJiaTAY2qxfe3HUXGTkl8ueT6+SNy7a2ubPuWjjURn76M9ye2Zhq+7j5eXloBHZprCiVs5fly7T88vl/HXpsrCi1vK7woraRts3dM1ZtOXfSFPP19D1hubQWTh7LgorauXLvxyUU+dvkS//fKBN3+fmAqRKO9Zvlz2xYCzgNcr0/WJgHfCodQPTg1aYfnQ3/XN+cZRm4OHhgYeHc6tLtjel1Vpm/XcXKjcFb1zdt81tHdcNjmHhhgye/mE/P907os1OSy3FvIO3rqVha1fv6jaU1mCPvaCh529MzXQucibgsmnOBfsKuHbN+3BpylJs+tqg9VoIoRRC7AHygd+llNsaaTdTCJEqhEgtKHDNuvDnIrkl1dy4aBsZBRXMvXYAkf7Od7uuj5fKjZeu6MOhk+XMX5fR5v23FLM6py1r0bsaHXPQPM6V+WrXE4sQYg3QUD56u9PCSCn1QH8hRADwvRCit5RyfwPtWlyPpYOG+e3ASeZ8uxe9QTL/xoFc0C203cYyLimcy/pF8d7aIyRF+THGQfVfnMk5oSt3Mh1z0DzOlflqV8EipRzb2DUhxClz7XohRCTGE4mte5UIIdYBF2PMWdaBk9DoDLz+6yE+2ZRJn2h/5l03wGkR9s3h1cl9yDpdyR1f7OT5y3tx/WDHVK10Fv9mFZe9dMxB8zhX5suVVWErOFPXZTrG4mJ1EEKEmk4qCCHUGBNiHmqzEf4/JOt0JVMX/MUnmzK5eVgc3951nksIFQAfDze+vH0IwxNDePL7/dzy+Q6OnGo4+6wr4+r1zM8lOuaydbR0/lzZeP8asEwIcRtwDJgCIISIAhZJKS8BIoHFQgglRiG5TEr5c3sN+N+IwSApr9FRXKXhl30neH/tEVRuCv5zQzIT+kS29/DOwl/tzqc3D+LzLVm8+3sa497dQL9O/ozvFcHQLsEkdw5w6VMMnBtxCucKHXPZOlo6fy4rWKSUhcCYBn6fh7FiJFLKvcCANh7av57bF6eSUVBBSZWG0mptnWp4E3pH8NykXoT7eTZ+g3ZGqRDcNiKeKwdEsyw1h1X7T/Lm6sOE+KjY8WSj2leX4VzRo58LdMxl62jp/Amjx+7/L0JCQmRcXFx7D8Ml2LVrF2p123tyuSJqtZqOz4WRXbt2kZyc3N7DcAmysrI6Phcmdu7cKaWUTZpQ/l8KlpSUFJmamtrew3AJUlJS6JgLIx1zcQZvb28qKyvbexguQcfn4gxCiJ1SypSm2rmy8R4AIcTFQojDQoh0IcRjDVy/XAixVwixxxSnMqI9xtlBBx100IERlxYsJqP8h8AEIAm4TghRv8jHWqCflLI/cCuwqG1H2UEHjiGnqIql24+xdPsxTpbWNP2C/wd0eHWdzbkwJy5rvDcxGEiXUh4FEEJ8jTHVy0FzAyllhVV7b1wspUsHHTRFjVbPKyv/4b/bjqE3eUqo3BQ8PTGJm4bGtvPo2pcOr66zORfmxNUFSzRgnbvgODCkfiMhxJXAqxjTvlza0I2EEDOBmQCdOzunDnsHHTSX/PIabv18BwfyyrhpaCy3DI9HbzDw0i//8PQP+1EIuGHI/z/h0lhxsOa+3tVzarWEhjy1bD1ve8yFS6vCMGYvrs9ZJxIp5fdSyh7AFcCLDd1ISrlQSpkipUwJDW2/1CMdNE11dTU1Nf9+VVBxpYabFm0nI99YauCFy3sTH+JNYpgvH09LYVT3UJ5fcZCDeWXtPdQ2x7wrX3PwFHdckNDsBfFcyanVEszR99ZzYut522MuXP3Echyw3qp0AvIaayyl3CCESBBChEgpTzt9dC1Er9ejUChcPlDP2UgpSUtLY9HXP1BSdObt8vBQY5AGNLVG4RIaFsmt10yia9eu7TVUh6M3SO5ZuovMwko+v3kQw+rVknFXKnh7Sj8uem8Dz604wDd3DP1Xf16KKjUs3pKFubKidfXFBeszmr3bPhfjV8xzUK3Ro1YpmD4s3u5ntvW87TEXri5YdgBdhRDxQC5wLXC9dQMhRCKQIaWUQohkQAUUtvlIG0Gv17Pk+9VkpB3idEE+er0OhUJZZ5GQUlp+ltKAr58/id2TuPqikYSEtF/xKkfw3Z87+Dt1G/mnTqDXas+6Hh0bz9hLJhES1lAuUuPcFJw6weL//cKxzAxGXHgRM6Zees4vsvPXZ7A5vZDXr+pzllAxE+zjweyx3Xj6h/38cSj/nEisaQ/1VTP1K2R6qdwsRbvMRdGgefaEcyWnljWLt2RaqofCmXlojPrz6EplF1xasEgpdUKIe4DVgBL4VEp5QAhxp+n6fOAqYJoQQgtUA9fINgrOycvL4535n1FeVsoZrZ2xayEUSGkAoG/yIEaNv4Tg0DAUiqa1j2UlJaQfPshLb39AeVkpw0aNZdoV43F3d3fSkzie/3z1A1vWrSEuoSuDhl9AeEQUbi0YvxCCsIgorrxuOjqdjq0b/uCWWQ/gHxjIg3fcTGzsuWd/OF5cxdy1R7ikTwRTm9hFXjsohvnrjOVr/y2Cxdr4PCUlxiJURiSGMDA2wLKzLqrUUFipYXhCsMXW0pi9wPr35j7OBfuK+ZRyoqSaVQdOAJDcOYCRXUOaPGG4shHfpQULgJRyJbCy3u/mW33/OvB6W46puLiYp195Gzd3dy6fcj1BIY612fgFBJA8ZBjJQ4ah1WrZvnk9t9/7IDGd43nm4XtQqVz3jyU3N5enX3mbHn36MfuJFxx6snBzc2PEheMZceF4igtPM/+LbynIP4GU0C2pN3PumGaX4G5v3v4tDQE8eWlSk/PjrlQwfVgsr6w8xIG8UnpF+bfNIJ2ItWpm8ZZM/jxcwPCEYN6/bsBZdgNzPfjFmzOJDvSiSqNj7tp0qjS6OuWezYtslUbP3uMldcoVtyfWTghrDp6yfDWPe75VzXszanelXWowV1b3ubxgcTVWrE/lmyWfcOus+wkJc/4O0t3dneGjxjJ81FjSD//DXQ8+RmBwCC8+9oDLpWL55rfN/PztUm6e9QABgUFO7SswOISrb7rV8vPfO7dx7c0zuWnmvVw2op9T+24NeSXVrPg7j1uGxREdYN/7d01KZ97+LY1vduTwwuXnvmCpi1GwpsQFnXUCqdLoiQ7wJLekhrWH8sktqWHa0FhGdw+lWmtg7tpDFgFjPtFUaXT8ebiA0d1DW2yfcSRmgbd0+zGyCqvYkFbA5oxCth4t5O2p/dmRadTae7kLArw8CPfzYHNGIctTcyxVRhs7fbmyuq9DsDQDg8HA0s8X8tDTL7eLWiqxe0/ue/QZjmdnceus2Vx1/XSuHje8zcfRGKtXfMddDz2Jh2fbJ6jsN3AISX0GMPeVZxjc9QXCw11TbbTkr2yklNw8PM7u1/h7uTM2KZyf/s7jqUuTULm5/qnMFtYqnOnD4vBSKS277oyCCl76+SAxgWqWbD1Gcmd/cktq0OmNKuZDJ8vYnlWMVm9g5sgu7MwuZlO6cXG+44IEMgoq2Hu8lKcmJrHm4CnLKcbcR1sJGPNz3HthV0Z3D7WcoBLCfBAC/jxcwOItmZgeiyqtZIZpDs7vFmqZj8Vbspi79ghVGj0PjOvm1DE70i25Q7A0gydeeotLr5za7raOTrFxPPLca7z1wpNcMnwAXl5e7ToegNOnT6P28nKYUDmdf4r3f9hIVVE+QqHg7skXEpvQ1abqyF2lYuzEK/j2t03cfdNVDhmHI9HpDSxLzWF8UgSdApv3nk0eEM0ve0+wPq2AcUmuJTSbuyBZq3DMu25zNLl5Rx8TaDzNJYb5UlqtI6PAmLfsaEEFvp5KNqUXIiVszihkRGIwVRodRZUa1hw8xZ+HCxja5ZSln8KKWtPirOOBcd2dNAvGeZi/PoO/c0o4VljFibIaqjR6uof74q92o7Rax9GCCpIi/U3CUOBp2iR4KGFDWj67c0p5fEIPq3mU9b6e6cvRdiRH2mw6BIudPP/2h4SGR9BnQJP51yzo9Xpe+OwnThzaia62BpphbzDodXTuP4IXZkxu8LpCoeCmGbN4/s15vP7so3bf1xmUlpYy68HHuOPBx5tsK6XkuUXfcyptNzpNbcNtDAbU/kH4R8QS1LkrBr2e97/5lcLs94gfPJbnbm94TgBqa2pQebhmSv+tR4soqtRwxYDoZr/2/G6hBHur+H73cZcTLM1dkIK8VSb7itG1tkarY8OR02QVVjEkPojhCcFEB6rJST3OkVPlZBRU4uuppLxGz+nKM56FCWE+pMQFsvVokcWbavqweKCu0Hr39zTTK5zrSWhtEzJzuqKWbZlFAAR6ubMpvZCBsUE8PqEHY5PCqdbqScsvp6RKx+6cUkZ3D61jMzE+jwAkRZUaixBxhuHekTabDsFiB0+/9h6+fv6MvWRSk20NBgPPLPyW3P3bAAjt0oukMVejUvs0q08pJUc2/cK9z77JvOfnNNgmMjqGnOzMZt3X0dTW1jLj3oe486En8A8ItNn26f98w7Hd64nsmULPMVNw97R/1x4Uk0jCsAkc/O1rnvqohpdmXd9gu1N5udx8dYPJF9qdX/bl4a1SMqp785093JUKJvaNZOmOHMprtPh6uo6HYHMWJLMX1Najp9mWWVznWpC3O9syixjdPZRarR6ASo0OALWbknKMv1MABsDTXUlqVrFl4QZx1gloSkrMWeo2R2F9Quke7kuNVo9KKdDoz5wuMgvOZIg+v2sI+/PKuKBbKMmxgSwwGe5nnt+Fg3llJEX6cueoxDonkCBvFV4qJa+uOsTe46W8PbW/RTCDYw33jrTZdKTNb4J3F/0XvU7HuEsvt9muqPA0j778LtraaiK69adTryEo3Fovtw+sWUZkj2SeuvGiBq9vXreGgvxTPPvQ3S26f2tSgldUVDDj3oe44fa7iezU+Adcp9Uy65Fn8A2Npst5F7XaU2zb0vdYNO8d3Bvwjvt75zYOH9jH6083LIxt4cz06Dq9gcGvrGVEYgjvX9ey2nQ7s4u56j9beHtKP64a2MnBI6xLS9Pm13f7XbwlExBMHxYHUCdeBSA6wJOU2AB2ZJWgVinILa6hRmdA7a6gWmtAqRCW/GnW+Kvd8PN0J6e4GoAh8YEM7RLCpP5RrDl4isIKDQs3HrXMd0PqInvVSQ19LurH3tjCX+3GqG6hrEsroLRaR3SAJ3HB3jw0vjs7sopsuk7Xj/MZ3T3UIlzaI1WLvWnzO04sNli5ZS/phw4yc7btRWrOW4soOp5OrzFT8PQNcOgYeo6aTOp386ERwTJ81Fi+/e/nLPjqB+64/gqH9m2LT/63inWrf2HanfcRFhHVaLuy0hLufvBRksZdg3+EY3K0dR5wPs/O/4ZX7rvprGv9Bg4h/0Qer8xbxBP33u6Q/hyBWQ12ad+Wl3NO7hxAp0A1P/6d53TB0lKsVTRm12AAL5XS4rE1JD4IkGzLLCa3pIZwvxryrLI5KwRUaw0IaFCoeKsUlFbrKK3W4e2hIMjLg+7hfsxde8TiahwbZDwNb0o/bfGwsjXW5uzUMwoquPnT7eQUV+OuAK3BdvvSah2/HTxFtalhbkkNuSU1iN/TGBgbyOItmUzqH21xQzYb7AsrNQSbTidvT+1vES7m53HlOJZz273EyezYsoFLJ19js80rS/+g9GQ2A6+Y4XChAqBwc8Nd7U1ZaWmjba6+4WbWr/nV4X03RmlpKetW/8L9T75oU6iUl5Yy64E5JE++02FCBSC8Wz/y0/c2en3cxCupKC9n3uLlDuuztZjVYBd0a3nMkxCCSf2i2Jx+mtMVDdun2puxSeGM7h7KoLggdmYbVV0jEs3BfsaT6qmyGh69uCfThsYSE6jm8Km6udDMsqQxXYrOStVUWWsgp7iaQyfLmD2mKzGBauKCvcguqmJwXCBD4oPILa7i3d/TLGnmzWqysUnhPD6hR7PUSUWVGm79bIflpNSUUDFTrTWgNB3Uw3xUxASq6RLixdy1R5i7Np1nf9xvlc/L+Hx/55Tw6qpDLN6SRZC3iren9q8z3ikpMc0ef1vhEMEihBgnhPhYCNHf9PNMR9zXdK+mCn0JIcT7put7TWldHELe8Rwio23vDLN3rafnaOd6IMX2H8lzH31pu02XBDIz28be8uTLbzPtjvuabHfPI08wcPKdeHj7OrR/hUKJXqfFYGj8r/qqG25my7o1uIKqV6c38Ov+k4xNCsfTXdmqe03qH4XeIFm574SDRudYzF5Z8/44wqb0QkZj6rueAAAgAElEQVR3D7Wooib1jyLI252swire/u0wmzNOk1NcTWVt06uztfK01kqwhPkYVUCR/mq+2JrFkq3GeJHR3UPp3zmQbZlFLNl6jLlrj/DQsj0W9dGrqw7x0s8Hm61Gmr8uneyiKrvaKutpfM3Dzq/QkFNcjae7GyMSgwEoq9YSF+zFoLggpg+L5/EJPegebrTLVpvsTPWTTzaUjNJVcNSJZRYwB7hRCHEh0N8RN7Wz0NcEoKvp30zgP47o24hsMpJbW1Pl8IWzPgHR8ZSePGazzejxlzB/yTKnjgOMzgmVFeVNZht49uPvCEvsi4ePcwL6InsO5On/fG2zTVhkFFVV9i0CzuSvo4UUV2m5pE/L1WBmekT40T3clxV7Gs3F2q6Yd9FPTUzi8Qk9LPYAMAqdokotAV5ulNcaXYijAzyJ9Gvai8/Ho2GBXF6ro2+0Hyv35lFUqcXHQ8nM87vQt5M/NRq9pV1MoNqiRpqSEmOJLWluxt+tR+1PQxgTVNc5xV9d1/KwfGcOt4/owojEYPbmlpFVWMU7v6dZxhjo7QGAWtV+FouWFhVz1IgLpJQlwMNCiNeAQQ66b5OFvkw/LzHlB9sqhAgQQkRKKVu9pQuLiOJ4dhadYuMabePm4Ym2pqpZHk7NRSCa3HlnpB2iW1Jvp43BzKf/+5W+Awc32e7EwR30vXS608YR1Wswe35YRL2cpHUIi4gmKyuLXr16OW0c9rBy34lWq8GsmdQ/ijdXH+Z4cVWz42GcjbVnUcIFPhRVanj398OA4IJuoSSEepNRUElJlVG1GxfszeaMphfr8lp9g7+v1hrYm3tGlealUrL2n1NkFFRa7CwAob4eTE6OZmxSOMtTc3hqYlKdWBd7KKrU2H1aAcgqPNNWAYT7euCtciOvtAZPN0FxlZbnfjpAuEmwxgSqSYr0bTR41HocbWW0b6kdx1Enll/M30gpHwOWOOi+DRX6qh8EYE8bhBAzhRCpQojUgoKmPTkA7rxpCpvW/W6zjZd/CNVlRTbbtJbSUzn4hdmOfTiRe5yxg50vWA7u3U2/lLNqrZ2FV1AYZfnHnTYOhULZZFyQu8q93eu6aB2oBjMzqZ/RrvXT366nDqu/wzUaotOZu/YIc779m4yCSuKCjQt+uK8HWYWV9I2278TvprAdieLroaRaqyejoBK1u4LsoipLX7uOlQDCEo3fkjovi7dkUVqta1E0jAFIy68kr7SGIfFBTOoXTVywF4PjAi3u0hN6R3DnqESL3aQxVVdb1ldpqR2nRYJFCBEnhHhTCPGdEGIRECOEsKSZlVLOa8l9G+qqgd/V37rbWwys2YW+oqOjOZ1/ymYbtX8wlcX2CaqWcuLQLuZMs+3xpVJ5oNE4vwZ2j9792L+7aZfc1x6+g7QNPzptHFJKpA0bC8A/e/cwYEDLXHsdxeb00xRXaZnYt3Enh+YSE+TFgM4B/Lgn12H3bAkNqUmsF72iSg1/ZRjr7MQEqskoqGREYgiD44LwV7txqryW3JIaiqt0lkh7W+gMjRv0kzsHEOzjQXmN3uKqHOjlzuD4IJI7G51qth4tbJHB/gzS6n/beLkblyWlOPO9mb9zilm28ziX948mr9TKCcMON3xzDrXZYxLrZIFuibrKHlpqx2npieVH4BBG+8c4oB+wQQjxoRDCo4X3bAh7Cn01qxhYc1EobO8yn7j5Mk6l/e2o7hqkovAE0Z1tp4fvm5zC0p9sn64cwW1XT7BLsHiqvYjsmcKx3RudMo6CowcIjrOdnkMoFO2e7fjnvSfw9XDj/G6Oratzeb8oDp0sJ+1UuUPv2xwa2jmbd7hjk8J5aNketmcZPcMm9I5k9piuSClZtvM4pdU6y2tCfVWM7h4GnDHGN0Sg2r3BXeSQ+CBSYgPJKqxCAAM7BxLo5U5xlZZlqcfxUhn/hrdlFrW4IiUYo+Aj/Oxb3syxWnppzANmxtNdQY1OonZXUFxZy6b000T5G1Vhe44VM39dhsUTrCHMrsjmYFBwzWqZLf2rU0opP5FSrgWKpJQzgAQgC1joqMFhVehLCKHCWOhrRb02KzDWYxFCiKFAqSPsK2ewvT/x9fNHp62lvMB5xlQhmn6bOscnkHss22ljMKNUNq2CMvPq7OmcOJSKXuv4ndSxXet5+V7n2XAcQa1Oz+oDJxnXKxwPN8eowcxc2jcKhaBdjfgNqUnMO1yzd9jwhGBmj+nKnaMS8FIp2ZxRSCerrM7hPh7sOlbKL3vziAlUE2JDsBRXa8/6a+wb7YeU0pKIUgKbMozOEnHBXsw8vwvPX96b2WMSmT2ma6tcc4O8VXQOss+mVak5417s62l874cnBDP/hoEEebtTrTXw017jMhVhEizbs4o5eMIcVtDYunN27jBXdDtuqWBZYyrABaYnlFLqpJRvAuc5ZGSmewLmQl//AMvMhb7Mxb4w1mo5CqQDH2P0UHMYgUEhTarD5r/5Igf/+JacvX85smsADAY9wo5dt3GH1DautTU11dRU22fE7DlmCgd+W+rQ/ktPZOMVGNpg5L01NVXNjxx3JH8eKqC8RsdlDlSDmQn19WB4Yggr/s5rN5dqW2oS82I37/pkS1be3OJq4oK9ePTi7pYYEzc34+pbWKUlp7ianJIzNjHr7Yuqvu8u4OEm0OmlaUGuGwtjDCb14ppBMQR6qfBSuTF9WNxZEe7NVSG9elVfogPsy0Vndi+OD/EhuXMAWr0BP7U7E/tEWYz3CaHepMQaUyENiQ+iS6gPQ+KDKK7U8u7vh8koqKgzRrMrsjknGpz9PjhTNWYvLRUsDwL+QohUIMpkGL9RCPEhDi4LLKVcKaXsJqVMkFK+bPrdfHOxL2nkbtP1PlJKh+bkmD3jJhYvmGczZkLl4cHi/8xFW1PJ1q/nUlFkWxA1h9ITx/ALazrKOvtoOgGBwQ7r1xavPv0on330rl1tH5s6Cq+AUDJ3rHVI33qthoNrlvHuMw/bbLfpj9/oP9hhe5wW8dX2Y0T4eTKyq3PKS0/qF8Wxoir25JQ45f6tof5itzw1hyVbs8kqrOK9tUfYbsrxJSW4WwmN8pozKjJrcanVS4ugCfNVoRRQq5NknDZuHsyVBEK9Vcwe0xW1u5JN6ad59scDPLRsT4OqopaokBJCffjp3pFMTbE/88He46XsOlbC9qxi5nz7N0u2ZlOjkwR5u/Pm1f1Qq9yYPSaRfp38WfJXtin2JtsUOHmgjmrMHpuHK6jGWiRYpJQG0yJ/PsbYkQhgILAfY1zJv4bIyEguu/pa3n7xKTLT02y2feuRO1jwzqukbfqZA2uWYdDrbLZvCoNBz6F13/P83WenLrFGSsmyLz7lxccfaFV/9hIdHU33Xn3ZuHa1Xe3fefI+pDSw6/uF1FY0nkGgKapLC9m29D36TLjR5mnleHYm2zev55E7209Vll1YyYa0Aq4dHIOb0jl2not6R6ByU7Dib9eMaTFTVKmhsKKW5M4BJHf2p2uYD74ebvh6KMkrrSHUp67dom+0r0U9BObcvmcETa3OgF4ahUmtzrjhU5s87jxVSqYPi+Oh8d1JCPWmS4iXJcfWlJQYk/tzGu/+frjFhvwgbxVvXN2PvtF+dr9GKSApwgd/tTtTUzqR3DmAxFAfftida7GZmM9oHkrBxb3CGRIfRFm10bhvDpK0B1dQjbUqjkVKWYXRxlHf7vGvYsr4EUwY1p83P/qU75YuIXnIMEZeOB63BpJMqr28WfTOK7z05Wq2/Pdt+l58g10nDmsqik5xfN9fFB3PoOeoK/ELsJ0qZsW3Szl/zEVtWrJ4zh3TuOH2WYwc03AOs/q89chdFBbk88Tr76OtqcLTNwCf4EjcPdWAQCgUePoG4B0UUSfgVEpJSe5RsnetQ0rJgvffxtuncffUdat/4fDBfXz+n7mtfcRW8dX2YygVgmsHOS6VTX38PN25sHsYP/19gscn9HTZAmDLU3NYuNGYFWJ091B+PXDmRJ8Q6s3whBCWbD1jH+wfE8SGI0ZPS3+1G+G+HqTlG08mvp5uFsO/SaYQ5e9JXmkNQd7GpJTmnXpGQSUT+0bVcd9dsD7DtJCDl8qtVTm23r12ALd9vqNOvAqAm0Lg7+lGYdWZFP96CSfKaimuqqCoUkO4nyfbMos4bHG+kKhNTga1esn2rCKKrEoENCdI0hUqS3YkobQTHx8fnn/kPqSULP5+NW889xgz75vTaHnip268CM2UUdzz5Etoa6oJieuByssHbXUlNeUl1FSUNFyPREq8AkOYc8vVxMY3/eH4/Zcf0et0zLzOdvZlZ5DUN5l/9u2hZx/7Ei0Eh4ax4K2XAGNyypO5x6mqqgBprF3z9cb95B3Yjqa6rm3EP6Iz7zz/OD6+tneIf61fS3lZKQvfe71lD+Qgymu0fL09h3E9w+vsvJ3BNYNj+PXASVbuO9GiOi/OxLree5VGBwgm9Y8i1DeD1QdOclFSBI9d0tPYWMCPe3KNxbBOV5BVWIWnm7AkmzS7EJtVZRF+HpwsM/79xAR5MX1YHIPignj913/441A+j17co45AMTMlJYYqjR6Qrd7RJ4T68N2s4WeVAdAZJJ4qJZgEi5dK0CnQiyxTCv2swirCTd5lpdU6EkK9mdQ/mkAvFcWVGv48nM95XYIJ8DKOW61SWLJDnyt0CJZmIoTg5skXM+Wikdx+zwPcNOMeomMa3pWqPDxY+NaLaDUa0g//Q2VlBT4+vgQEBeEfGIRa3fKoaSkl3yxZhJ9/AM8/0nTeLmfwyF3TuXHG3ezYsoELJ0yiU+c4u1/r5x+An3/dk9jAoS0vs6zTatmyfi1ffeLAjD4t5PPNWZRWa7l7dKLT+7qgayiJYT4s2nSUy/tHtbokgSOxjtq2rtyYEOpDabWOrMJK7vpyJ/06+XP/2G5MHxZnKefrrjxiSUmvVAiqtQb81W5c3i+aQG93th4t5GRZLQFebrwyuQ8JoT68+3uaZXGf98cRPrtlcJ26LEHeKoK8VQ4t8Wu+X1FlHGPfWWc5ZYzpYdxw/vi3UVimnTLbggTJnQN49OKerE8rYOORAnYdK2HxliwCvVQEerszoXckCzceZfaYrk4vR+wsOgRLC/H29ubz+fO464FHOX/sxQwYNLTRtu4qFT379HNY3wWnTvLloo8YOeYibpvSfkWt3N3d+ebzhRQXF/PMa+/h7evL5VNvROXhyFCmppFSsmThPCZf3/7ux6XVWj7eeJSxPcPp08k5edKsUSgEM0bG8+j/9rHmn3yXqi7ZWDGqKSkxbD1aaBEc2zKLCDbZWcy/e2piElWavWzLLLakzi+t1hEdqGZKSgzVGj3uSiXPX96LhFBzET1ju9ggL56aaEwp2Fap5YO8VSyaNogHl+3h/K4hTB9uFJLW8TpgPM1szypmR1YRD4zrxs5sY9T9+sMFlnQx5sSU1Rp9HaF4LtEhWFqBh4cHn3z4Lk+8/BaHD+7j2ukznNpfeVkp//tqMXqdjvffeAk/P/uNh84kMDCQea8/z/I1f7HwvdcJi4xi8nXTcXN3fpXD4sLTLFk4j+EXjOXqMY0L97bio3XplNXouH9s1zbr86rkTizYcJTXVv3D6O6hTnMWaC6N6frNKeDnr0vn7+Ol9OvkbxE+ZoHTNTyHk6ao9JhANaN7hBHo5c6UlBiLzWb2mETWHDxFYIrxJDJ9WDxeKrc6C7EzKi02RnJsIOvmjAZgwfoM/jxcwIjEEMqqNezNLWPqwE4EeKs4kFvKoLggFqzP4MFx3XFXHuHeC7uyPq0AkEzqH82KPXmkZhVZ8qi1t82kuXQIllYihODVp+Zw/+PPUVxUSGCQY1x+DQYDtTXVFBcVcmj/Xg7u24ObmztPPDCLmBjXCYSyZsrY85gy9jz+t3YbH775IiFhEVw06apG7VAtxWAwcOxoOr//8gNCoeDdl58lOLhtXK1tkVFQwaebMpkysBO9o51/WjHjplTw2MU9mPnFThZsONomKrjm0FBFxOWpOWeV4QV4e2p/lqfmsCHNuINPCPVm+Z3DzrKTAFRp9Ly66hBbjxZasijXX4Dby5BdX6CZn395ag6bMwo5acqb9viEHnx2izGpa3LsmdLe5mDSuGAvCitq69S7PxfoECwOYs49M5jz9ItEREUzeNj5dO/Vx650IseyjpJ7LJsTuTnkHT9meo1ACIGnWo1/QCA9evXlvluvM0a9nwNcNWYIV40ZwokTJ3jjo08oKSwkIjqGnn36ERkdg39gULOfpaKsjI1/rOZo2iEUSgWd4xN5/dnH8PdvuwXcFlJKXvjpIJ5uSh65uEeb9z8uKZxL+0by7u9pjOwaQt9Oji8611Lqq6PMP1sLBDNmQTA2KZxnf9xPUuSZ99daQJnr2psrRj60bE+de7VH2V5r6gs08/fWasARiSEUVhizP08fFl9H6I5NCre0W7gxkwN5Zcy7PvmcES4uK1iEEEHAN0AcxlQxU6WUxfXaeAIbAA+Mz/KtlPLZth2pkejoaL76dD5FRUV8uuwn/lj9M1JKvL19CAmLQG/QExvXhYwjh8k/eQKQaGpriUvsRmx8ArddewUxMTEuZXxtLZGRkbz74lMAZGdn8/3av9i/ZydlJcUY9Hr0BgNSSoKCQ+jRux/xXbsTEBiEpraW48cyOZWXS2Z6GkWnC/Dx9WPEmPE8ed8Ml5yj5anHWZ9WwDMTkwj1bVsbExhPzi9f0Zvd2cXMWJLKj3ePcLpHmr3U371bL66NlQ1OCPVhZNdQXl11iGAfVaMCqaGSvdB2tpXmYh7z8tScOqWbwejgYD3ut6f2576lu9mUfprNGYWNzpUrIlyhwl5DCCHewJiH7DVT5chAKeWj9doIwFtKWSGEcAc2AbOllFtt3TslJUWmpjo0QL9RysvLOXHiBEqlktVbdjPxgsEuJUBSUlJoq7lojPz8fJb/tonMI4cpLSnCXeVBp85xXJDck27dumFvNurW0tK5yDxdycT3N9Knkz9f3T4UhaL93tt/TpRx9X+2EB/qzdIZQ/H1bJmdy9vbm8pK56XEsedEYd0GYPGWTHZml7Ap/TSPT+hhWWQbupcjTyzO+hspqtRYBIfZA6whtaEx6t5oe1lz8FS7GvOFEDullClNtnNhwXIYGCWlPCGEiATWSSkbTWcrhPDCKFjuklJus3XvthQsro4rCBZXoSVzUVSpYfJHmymt1vLTvSNcovDWn4fymbEklT6d/Fly6+AWCRdnC5bmsmC9Mevv7DGJmGPxzeojZ+PMv5HmCFjzCcdaqLY19goW13AfaZhwc5Zi09ewhhoJIZRCiD1APvB7U0Klgw4cxdGCCq6ev4W80hoWTU9xCaECMLpHGB9cn8y+46Xc+Ml28svbt9iZIzCnKTF6fimZuzbdpdLEt5Tm5P4C0e6pWuylXW0sQog1GPOM1edJe+8hpdQD/YUQAcD3QojeUsr9DfQ1E2NeMzp3dl6ajQ7+/Wj1Bv638zivrPwHN6WCL24dzMDYoPYeVh0u7h3BRzckc9/Xu7ls3ibemtKPkV3bRqXoDKyN4W3pQuwKWD9vh/HeDqSUYxu7JoQ4Za5db1KF5TdxrxIhxDrgYozJMOtfX4ipVkxKSopr6v86cGmOF1exat9JPtucSV5pDSmxgbx7TX9i7KzR0daM7xXB97OGM+u/u7jpk+2c3y2U6efFMqJriMPrw7QlrpALqy05F5/XZb3CMCa2nA68Zvp6Vp1bIUQooDUJFTUwFmjfRFEd/OswGCSXfbCJA3nGmh9D4oN4+co+jOoe6jJOGI3RM9KPVbNH8unmTD7bnMVti1PxdFcwKC6Ivp38uW1El3NmF9zBuYMrC5bXgGVCiNuAY8AUACFEFLBISnkJEAksFkIoMdqLlkkpf26vAXfw70ShEFzQLZQrB0QzukeYVQqRcwNPdyWzRiUyY2QXNqQVsPHIabYeLWTB+qPMHHlu7YQ7ODdwWa8wM0KIi4G5gBKjQHmt3vXLgRcBA6AD7pdSbrJ1z5CQEBkXF+ecAZ9jZGVl0TEXRnbt2oVarW664f8D1Gp1x+fCxK5du0hOTm7vYbgEO3fulFLKJp2+XFqwmE4iacA44DiwA7hOSnnQqo0PUCmllEKIvhhPLTZDnzvcjc/Q4W58ho65OEPHXJzB1Vyv25N/g7sxwGAgXUp5VEqpAb4G6hQekVJWyDPS0Zu2KvzeQQcddNBBg7iyjQUgGrB2Vj8ODKnfSAhxJfAqxliXBvPId7gbd9DBuUONVs/89RksTz1OjVbPqO5hPDCuq8vECnVgG1c/sTTkcnPWiURK+b1J/XUFRnvL2S+ScqGUMkVKmdJWKUI66KCD5lNeo2X6p9t5b80Rukf4ckG3UH7Zl8fYd9b/K4Ii/z/g6ieW44B1FFQnIK+xxlLKDUKIBCFEiJTytNNH5wAqKirIz8+noqICKSU6nY6oqCgiIiJc3pXVkdTU1FBcXExZWRmlpaXsSMtBqVDSJzYUNzc3evXqhY/PueWNZQ8ajYbn3nyfotNnf1yllEgkCqGgtqaa/oOGcs/0qedMluuWYDBI7l26m53Zxcy9tj+X9zeWW37oou7MWf43c77dS3ZhFQ9f1Gh2p3MaV8sN1lJcXbDsALoKIeKBXOBa4HrrBkKIRCDDZLxPBlRAYZuP1Iqamho++99K9uzYSsOHLmM2WiklarWakLAI1F7GI76bmxunfv2T3GPZDBh8HrNvve5fJ2Cqq6v57o+t7Ni8gYrycoRC4OHhia9/AD6+fsZ/fn4gJbuPnkCn0/HfH1ZRXVWJlJK4hG7MuXM6Hm1cqbK5/Lg+lV9XfIdeb6wiKKw+C9J08K4sL2fqtNuI7WK7hoqUktS/NjH9jnvp2rMXD8y4yWUKvTmSxX9lse5wAS9d0dsiVACiA9QsuXUwT/2wnw/+TCcxzIcrBkQ3fqNzlMVbMi0Zj3/ck0dWYRVVGn2dEsXtXRLAHlxasEgpdUKIe4DVGN2NP5VSHhBC3Gm6Ph+4CpgmhNAC1cA1sg1d3aSUfL1yPds2r6eyohyFQom7uzv9Bw1lxr0P465q+Ru/ed0aps28h3lvvkxAgOvU17AXKSVVVVWUlZWxcvNuUv/ahFajQeXhQZduPbjiumln1b1vjPPHXGT5/p99e7j/ieepqqyg/6Ch3O9iwlej0TDn2VdQurkx7Y57UKtbbxcQQjBo2EgGDRvJ0SOHeeKlN6iqrCQuoRu3TJnossXfmkNJlYb31hxhZNcQbhhyth3UTang5Sv7kHaqnOd/OsDoHmH4q51fpbRtMX6OY4O8yCqsMv2u7nLmqiUBrHFpd2Nn4Sh341827eGrzxbQb+Bgzjv/QvycUHSqrKSEj+e9RVLf/jwy61a7ioc1B2e4lRYVFXHLnfcSEd0JD081Pr5+hEdGM3TkqFYJ2vqYd/HrfvuFK665iWsuGtGq+7V2LgwGA1/+tJZfvvuGm2bcQ6fYuFaNpymklBzLzGDzn2soLi7krReecljhs/ZwN/7wz3TeXH2YVbNH0jOy8dPY/txSJs7bxJyLurdJtcyWuhtbq7XqF/JqLMU/YCn0tWJPHg1lcXbWicWe+9rrbuzSJxZX5rUPFpGTlcn9jz+HyokqGb+AAB56+iX27trBjbfPYtS4Ccy4dpJL7dDr88WPv3PFtTfRL+UsBz6HYt7FDxw6nC8XfsCRQwd5avZMp/bZELW1tTz5ytsUnS6gz4AUHn72Vdzdnb+TFkIQ2yWR2C6JFJ0u4L5Hn2LE6PHMuOYyp/ftaPQGyVfbjnFel2CbQgWgd7Q/53cL5fMtWdxxfhfclK7pg7Q8NYe5a48AUFylJdBLxdajhWzLLLKot4oqNdz1ZSrbMosprKjliUuTLKcQa/WXNc7KHebIk5BrviMujE6n4477H8E/IJCZs+c4VahY0zd5EI889yrVVZXcNONuvlm1oU36bQn796QSFdN2Lt0KhYJpd96Ht48Pj730Vpv1C8ZTw/Q77+X8sRdz32PPMvqiS9tEqNQnKCSUB558kaNHDvHK+ws51zQRW48WkltSzQ1D7fvc3DCkMwXltWw84po+OkWVGgorNfh4GB0tvt99nLlrj7AtswiAao2eBeszeG9NGtsyjYVxD54od9pYFqzPoKhSY7OduTSBddZoe19bn44TSzMoLS3lrvsf4ZrptxOX0NWu1xgMBk7nn6K8rBStRoNWp0VbW8vy7Rnoamtw91AzfUw/OscnNKmLF0IwavwlnD/2YpZ+tpDszAwemXWLIx7NYXz6v1WERUQRGh7ZotcbDAZe+XI1Jw/vpqa8BIQgutcQnptxZZOntAvGTWDdbyt55PnXef2ZR9rkVPf1qg0MHDyMznFd7GpfXlbKI6/Oo6a82DQ+6zFKFG4qQrv04smbL8PHt/nG+Wumz2Drhj+58fZZXHXDzUy+0LmnRkfx24GTeLorGNMj3K72o7uHEejlzre7jjO6R4Olmtocs+qrWqPn7+MlFiECUF6jx1/tRmyQF/1jAjmQV8rmjEICTDYipYDbR8Q3WEGytWove08iDZ2EWnqK6RAsdvLjuh189dkC7n74SfwDAm22ra6q5JHXP6T89EkUSje8A4JRefuhdFOhUCpRuqtQqX3w9g9BW1vNopV/UXLia3S1NSjd3Inskcyzt0/Gza3ht0ehUHDDbXfy30/m89XPf3L9xNHOeORmUVxczIvvfIiUkhtn3N1ou6rKCua8PJfaihLO9pgz7rIDorsQP3gsXv7BGHQ6cvZuZtrtd9Fr3DU8dq3tZx01/hL++PVn5i1ezn03T23lUzXNkUMH6Zcy2K62z3z8HZk71tJ7/HX4BDdUhgh0mhoKjh7gweffQlNZhrvam+jeQ3j2lsvsFpRDzx/NwPNGsPyLTzl8cB+P33O73c/THkgp+f3gKUZ2DUWtss+VWuWm4PL+0Xy17RilVVr8vdrHiJ9RUMFLPx/k3gu78s7vh9mUXtchNeG0O28AACAASURBVCnCl/JaHRqdgVPltezNLSPYx4PNGYUkhHqTUVCJ2l1BtdbAok2ZAGxKP21RlTlCPdWa+jUtfW2HYLGDl96dT3HRaR5++uUmVV8PvvIBRccz6DHqcvzDm/9G6rUa8v7ZycwHHkMoFLz2xIOEhje8CF13y0zefO5xrhwztF2SJ0opWfjNT2zbtA4vL28umTyViKhOjbZ//Zv17P9tKX0vuanRhbU+Cjc3YpMvIKbvcPb8/BmPncjitQdsn9JGX3Qp8157HtpAsKQd3MdlV1/bZLtnP/6eE//sZMi199sUEG4qTyJ7DCSyx0AANNWV5O7fxs2z7sfD249XH72X4NCmd+ju7u5cf+sdzH31OaSULm2TO5BXRl5pTaM2hca4KrkTn2/J4pd9J7i+AS8yZ2I+SfyyN4+9uWVsOnIareGM+tFdKdDqJdnFlVTWGojwM64bfaP9iAnyYnhCMA+N7876tAKKqzRk5FcQ5e/Jsp3HTXcw3ssRRc1aY5Np6Ws7bCxNsGLDTsrLSpk28x6bQkVKyYwHHkftF8iQa+5pkVABULqriOl7HoOuvos+F13PIy++ydMLljfYVqFQMO2Oe5jzzMst6qs1lJSUcN2td1JTXc2dDz7Orfc8aFOoVJSXsf+3pQy97gG7hYo1Cjc3kq+YgU5Twx0PPYVWq220rRCCkLBwcnNzm91Pc5j3+TcMGTmqyUX75a/WkL17A/0m3tzsBV6l9iZ+0IUMnnov3UZO4vHX5zGziee3ZsiIUby1YEmz+mxrfjtwEoWAMT3tU4OZ6R3tR9cwH77bdbzpxg5m/rp0Xl11iIyCCoA6QgVAqzfFKdUaADhZVguAn9qdJX9lszmjkHl/HAEkS/7KRuWmIN10L39PNyaZYnjsKV3sinQIlib4fukSrrz2xibbzXriJaJ6DaZTn6EO69vD25dBV9/Fsd0bKSstbbBNZHQMai9vVm7Z67B+7eHeR57izgcfY+SY8Y2q7MxIKblnzlMMuPx2FE20bYqEIeOJHzSGW+64m6LTBY22G3PJJN792HkLqkajYdum9Yy8cLzNdi99sZrD639k4OQ7Wn1qUPsF0u/S6cSljGb6zFm8tqxpB46hI0fxz749rerX2fxxOJ/kzoHNXjyFEExO7kRqdjHZhW2XfTijoIJvTKllNDqjAHGvp8EL8z17ExoX7MXtI7owJD6IcF8P/jxcQHGlltHdQ/nzcAHuSuPno7RGx4o9uRRVanhl5T/c8PFWdmUXt8iI3l50CJYmMEgDai9vm21eWPwLbu4qIrr2dXj/Qgh6j7+Gx9/4oNE2V10/je+WLnZ4341RVlaGn3+A3cGND73yATF9h6H2s22bshf/iM4Mnnov9815Ap1O12Cbwwf20avvAIf01xDzPv+G8ZddabPNi1/8SvpfvzJo6j0olI7TOvtHxDL0ugc4vP5Hnvzo6ybb+/kHUNrIxqS9yS+rYX9uWYsN8FcMiEII+G6Xc0+n1jzwzR5Kq42fO/NJRauv2ya/vLbOzwLIKqxi8V9ZbMss4pTp+p7jJaTnV9A32o/u4X4MiQ8CoFpr4KFle1i44SibMwqZ8+3fvLrq0DmTK83lBYsQ4mIhxGEhRLoQ4rEGrgshxPum63tNaV0cRsrQEfzxq+2ilOl/rabn6MmO7LYOPsERVBTlN3pd7eVNSGg4P23c5bQxWKPRaPD28bWr7Sv//Z2q4nyikgY5dAzunl5E9Egm++iRBq8fPrCPKeNbFzBpiyP/HKBn736NXtfpdBz683tSrroLhcLxub2U7ipSrr6LkrwsHn5joc22fZMH8eWK3x0+BkewLs146hzdvWWCJdJfzfCEEL7fndsmLtZFlRoOnSxr1mvcMFpM1O4KnpqYxOwxXZl2XizDE4IpKKslp7iavbllLNmaTWywF49P6IHaXcGfhwuIDvAkyt+TgZ0DmT0msVW2lrbEIYJFCHGnEOJjIcS1QoifhRB3Oei+SuBDYAKQBFwnhEiq12wC0NX0bybwH0f0beb2ayaRldHw4gVQWVGOp48/wsER8fUJ65LEs4u+b/T6wKHDOLS/bdRhfn5+VFY07XNfUlzEoT+/o88lNzllHDVlJY2emrRajVOTVmq1GptZBF5aspLoXoOd/rlIGnM15QW5vLJ0baNtevcfyJZ1jV9vT9Ydzifcz4OekfZtVBriygHRHCuqIjW72IEjO5uiSg03LtpqUX/ZS2yoN0oBD43rxpqDp5jUP4roADXzrk9mcLzxFO9m+phszyzijgsSmD4sntHdQ8ktqSGvtMZk1BcsT805J9RhjjqfXwhcA2yUUo4QQsx30H0thb4AhBDmQl8HrdpcDiwx5QfbKoQIEEJESilPOGIABoPB5uJwIvc4vqFRjujKJrHJo9j5/UKgYfVLSVERg3rGOX0cACqVCq3W9odbSsnsR55m4OQ7nbJjB6gsPtVgvExtTQ3u7s41djZlLynJyyK0S/09kHPoeeFVHFyzDK4b0+B1lYcH0Z3jyM/PJyzMNWI+ALR6AxvTTnNp38hW2Z8u7h3BsysOsHT7MQbFBTlwhHVZnprToiDGjAKj/WfRpkxOltXy5dZscoqrASxfA9Qq1Col53cNoahSQ5C3iren9rfExahVCkC6fI4wM44SLIWm7MKvm36utdnafuwp9NVQm2igjmBpaaEvDw8PdDY8cBQKBUiD3fdrKUo3N5sCrriokOBe8U4fh728+PnPhCX2QeXlnFNDReFJvAMbXiTzT50gPMq5mW8VTaSuD4yKp/BYGgGRcU4dB0BlcT5qP9sLqkIhXM7lODWrmPJaHaO6t64+kreHG5OTo/l6Rw5PXZrkNA+qKSkxbM8sZO2hxp1GbHHK5BmWU1zNkPggqjQ6zG7Fpys1xHl6sWTrMTxVbgR7q5iSEnNWVmMvlRtjk8JZsD7DpbMbO+qcPhdASvmT6efvHHRfewp92VsMrEWFvpRKZaMGYoCQsAiqSpyfVqKmotTmwlBw6iTR0a6TRrwwO42wxD5Ou3/mjrW88MCMBq9FderMiePHnNY3gKamxub1Z2+9jPz0feh19rkFtxQpJYf++I7XHrnTZrvSkhKXy5C9ct8JPN0VjOza+sJ7Nw2NRaMz8M0O5xq3D51s3onFeoG1XpQ0Or0pPb4gJlDNgBh/Szbj3w6c5NVVh3ho2R4yCios3mBm1+M1B0+dZchvaeoVZ9EiwSKEiBNCvCmE+E4IsQgYK4SINV+XUq530PjsKfTVrGJgzUWv19u87ufvT3VZiaO6axCdppbU/83n3ReebLSNEKLJsbYlSncVeq2jDq51kVJSVVpIUEjDC5JSqUSnbXwz4Ah69OnHrm1bbLbpNe4a9q9e6tRxpG9ZRezAUU2mA9Jqatslh1lj6PQGVu0/wZge4Xh7tF5x0jXcl6Fdgvjvtmz0BucY8Zen5pBbYntDUR8DZ3a+AvBSGX86erqS5M7+7DpWQk5xNYPig5k9pivDE4LJKqwiIdSbPw8X8NLPB+sIkaJKDVUaPbPHJFpOLuZgTVfyGmvpieVH4BBGw/o4oB+wQQjxoRDCkVkZLYW+hBAqjIW+VtRrswJjPRYhhBgKlDrKvgKwf/9+Erv3tNnGOyiUwpx0R3VZh9qqCrYufY/+l93caFp+KSXHj2W51I70gRsncuKQc7zUju3ZSFTPxjN32zphOopHZ93apLfg49eMQqFUUpyX6ZQxaGuqKMk9ykt32c4woNPpmlTdtTXbMos4XaFhYt+W5ZRriGnnxXG8uJp1hxv3oGwNU1JiSO5s/BuztXDWn2lp9bVKI1G7Kyit1pF5usrSRu2u4IFx3Zh3fTIzz+9CiI8HM0fG89TEpDqJIc0Zk71UbnVOLg0lkGxPWrpVUEopPwEQQhRJKWcIIdyAB4CFwHRHDM7OQl8rgUuAdKAKcGhWxsXLV3D51Btstnn/+UeZPvNuki+/DbW/44yH5adPsOfnJXz0zmsEBDZ+31U/fsuocRPadEfalL4+JjaesvzjlBfk4hvqOBVdcV4mBRkH+OSDxrMYZxz+h65JvRzWZ0MIIYhP7MaJ3Bwioxv/Y/7g5SeZPvNuzrvhIYd7iO1f/RWvPzOnyXbZGUeItzNpaluxdPsx/DzdHJpAclxSOOF+Hiz5K7vZUfz2EOStYtH0QSxPzSHK35PZ3+yh/uFICTSmN0iK8CHQ24NrB8XwzIoDFFdpSe4cgLtSQXGllld++QeA3w6eJKuwCi+VkoRQHxIuOGOnbCjFi9nW4koG/ZZ+0teYFnwwCWQppU5K+SZwnkNGZr65lCullN2klAlSypdNv5tvEipII3ebrveRUjq0OlFFeZnNRR2MqpeF77/Fzh8XUZzb+t2pNBhI2/Qzhzf8xOfz59nsPzM9jcz0NGZed3mr+20O9sQMLHzvDdI2/sShdT9gcMAp4tiejWT89SsL575us93mP39n1o1Xt7q/pph+9US2b7Yd/e7m5kaPUVeye8UnDo2zyPhrNX4RnQmPbFpob/zjN2Ze77w4q+aSX1bDr/tPMiUlBs/6IeutwF2p4LrBnVmfVuC0SHzzAp5XWnOWUIHGhYqHUuCmFGzOKOTt39MorjLa3grKa9mWWcSSrdks3HiUhRuPWlRhT000ehU2Zj9x5XQvLRUsDwL+QohUIEoIMVMIcaMQ4kPaud68o0nq07/JxQPAy9uHxQs+JHv3Bvat/hqdpnm6WACDQU/Wrg1s+epd/CNi+WTu6zbzk21et4aVPyznP++81uy+Wos9i6S7SsWi998kLKE3O/73EXtXfUnpyWPNWmC1tdVk797A1qXvIg2G/2vvPMOjqLoA/N7dTS+kEBJSKKGE3jsoShEUP7CBiIhKUUQUu6AIIiKgIAIiSFUQUMSGNBEUpQqh917SG5BeNtn7/dhNCGm7SbYB8z5Pnmxm7txz5uTunLntHJbMnYHayLBORka6VfLB16lTh8irl42W+2BIb4KbdmTvqs+5HP432ekVy7uh0+URe+YQ//0wB7WDI5+PHWX0GiklyTeu4+vrWyGZlmD1vghydZLBHWoaL1xOnmpXA41K8N3eK2avuzD924Qwpnu9guCS+TgZnqhV3RzxcFbTLcwPZ40gO09yNEr/f7+clEH72j7U8nUtWG5cmFq+riwc0gZvV0e+/ucCC7ZfYOqm03y7+7LdzaWURoWGwqSUOmCKEGIW0ANoAXgDx4HSZ5hvQ8YMG8S7k6Zz/PABBj43oswd5xqNhq8/+4ip3//Fwd+WonFwpEaLLvjUqFfiXo7cnGzSkmK5EX2JxCtnydNmE9y0IysWzjM61PTjiqV4enmzcPZnlb7HilC7bhinjh2hYdPSd5/n8/4zveCZXlxLTGDy/JWc27kBoSqaiwQQAoo4HbWDIwFhrVg2f45JaZnXLF9C2073lOdWKowQAqnTkZ2VhZOzc5llJw1/BN3Qvnz0zXrO/PMr2qyb4+sIgUCgdnREpdYPZ+bmZKHLyy1kDwkI/EIbs2TuTJOHPbdu+I3WHTpX4O4sQ2ZOHiv2XuHe+n7Urlp2qKSK4O/pTK/GAawJj+SNnmEmh+EvLz5ujrzesz59WwTy6LxdpGTpe+SGmJMkGnoXO88lkKMDNycV6dk6HNUQ5u9B8+AqLNxxjS51qxLq58bfp+OJuJ5JiLcLl5MyWHc4CldHDVM3naaWb/7CDHnLUFhpuVoslbq4PFRqOYaUMgP95HnRCfU7BpVKxWeTxhEVFcXEaZ9Tu049+jz2ZJkP/nEDu8HAbqSnpTLxq1VcPvhP0eclABoHR9x8quEdWJvJLw82ORvludMnyMvLtWmSr7deHMJzL71KTnaWySmIfar6MeuD1yyiz7nTJ/l9zUru6dGbFwf2tYiMkniw3xNs/PVHHh1oPLqASqXiw6F9geL66XQ6srMy9avZBDg7u5S5s98Ujh0KJzYmipllrCa0Nqv3XSUxLZvRFsxV/0zHmmw4FsPvR6MZYOHJ7Dp+7gxsV4OF/14E9PMoHi6OXEpIIz4thxyDo3FzdCA9O5ucPDgWnUq3hgGM6V6P/Kn9iOuZ1PFzo3Odqizfe4UDV64zqV8T9l5M4u8zCdwf5seznWrfMpfy9T8XStwwac4UwxVFycdiIkFBQSyeO5NF369jzrRJDH35dTw8S16llY+buwcz3nnR7Lr8tmYVy+bPNnu95cHR0ZGVi+cz+YsF7PhrC17ePoQ1bkad+g1KXQZsCWIiI1i7chlBNWqybP5sqy+pfbxHR9atXUVmRrrRYKVloVKpKnV9PjqdjjMnjrHn37+QUjJn2qRK12kusrR5LPjnAh1CfWhX23I75NvX9qG+vztLd17iiVbBqFSW3Rg6smsdXBzU5DuJ2dvO0yy4CvFpObg6CDK0siBysbNGxYA2wfRtEcjH60/y95kExnSvVxDh+OFmgQWf1x2OplmwF82CqxQ4lcKUlqvFHDlcKoviWMrJiIF9eaR7R8ZP+QwhVDRu1oKadeoRFFLTaPh4c/DL9yto27GLVWQZQwjBhNf1YeESExP5aetuNv36I/GxMfTq+xiNLBhdWErJioVfotFomDt9skXjghnj6WGj+GrmVEJq1qJd567UqF3HpGE7c5GZkc7va1eTEBcLQhDWqCmT3n0Nb2/zRJM2F6v3XSU+NZvZAy3XLkDfLkd2rcMba46w9VQcDzQuf/6f8pA/LAY3d8f/dToOgKbB3nQI9SXmRib7Ll/jzZ71+SE8knE/HWXf5euGnkgtgILhq/zPGTm5zN52nvtLiUxQ2kowe1ghJqwREdTeaNOmjQwPr/zisaysLNZu2cmVi+e5dOEcTk7O9O3/VJkJrypKXGw0KxfPp1PX7rw46BGz1dumTRvMYYvC5OTk8PmiFfo5mGYt6NnnEbOGE0lLTWHhrOn06NOP5x7tZbZ6K2uLiIgIvv1pA1cvX0QlVAWLFO574CGT5qLKi06n49+tmwnfu4uBz47g4S4tzFa3udvFjYwc7puxnUbVPVk5vL3Fw8vk5uno/vk/eDhr+H10l0rJc3NzIz29fKvM8lMWj3+4UcF+EwAfNweupetXhDULqsL9DfyK9Uby50h6NPIv6NWMe7CBzZ0FgBDigJSy9E1kBmz/2nsb4+zszOC+PdCvX9DnfZ86+2tuXE9Co3GgYdPmODu7EFSjFgGBQUZXM5XG31s2cvbkcb76fBpubuaf8DQ3jo6OjH15GABLftrEjEnj6PPYk5XuwWSkp7F53U+cPXGc+bOm29VKJ4CQkBDGv3ZraJXs7Gy+WLyS9T99z4Ahw6gZWvm5BSklO7b9wb7dO+jaozcrFn5pd3HAijJzy1lSMrVM+F8jq+iqUat4+f66vLP2KBuPxdLHjBsxTaGOnzvLnm8HgHcbRzJy8gi/fI1dF5LQqAS5OsnV6+nM3pbM0chkZg5oUeBcCs+RzBzQ4paezO2C0mOxEBkZGRw9epTjVxOJunqZmKhIpNQhDG+yhb9ceXm5qFQqvH2q0qRFK5xdXMhIT+fGtSQunD1NFW8fPnp3jEX0tESPpSg6nY6Jn83l+rUknhnxcrkmpXU6HaeOHWb39m1otTk8/PhAHulq9IWpQljSFlqtlg+mzSIuJpr7HniIZq3aGn3AZmdlcen8WZIS47melEh0ZARabQ66vDxate/Ey89Ybq+OOW2x7VQcw74N57lOtfiwr2U3rhYmTyd5eO5OUjK1bH2ja4VXiFWkx1IS+T2RMH8PJm84yQd9GrF45yV2nk9kTPd6twyn2etqL1N7LHbrWIQQPsAPQC3gMjBASnm9SBln4F/ACX3va62UcqKxuq3hWCpCXFwcazZtJzsrC1c3N9o2qElQUBCBgZYLy28Nx5LPbzsOsnrJAnr3exwpJSeOHiI9Vb+2X61WUz24BjpdHqnJyWRnZ5GSfAONRkNYo6a8OOhRPDwqnrPDFKxhi9zcXL5YsoqjB/cT1qgJteuG4eTkRPXgEM6fOUVcTDQRly+SkZ6Go5MzdcMa0qlpXapWrUrNmjWttjjBXLaIupFJnzk7qF7FhV9GdTLrhkhT2HfpGgO+3sMr3ery5gNhFarDXI6lJGb9eYbZ284zpntdXu9Ztn75q8DyQ7fYwsncCUNhY4FtUspphsyRY4F3i5TJBrpJKdOEEA7ATiHEJinlXmsraw78/f155bknba2Gxeh3Tyv+13kB079airuHJ2NHjyiYYNZqtVy6dAkHBwe8vLxwdnbGzc3N7od4yotGo+GtF4cAQ/jhjx3EREaQnZXFP1s3UzesIb06tSR08OM2XYxgLrR5Ol5ZdRBtro6vnm5ldacC0K62D4+2DGL+9gv0ahxAk6CyV3Jam2c71cbVUWPSUFfh1V72sKS4LOzZsfQD7jN8/hbYThHHYkjulWb408HwY59dMAVAv6x23OjhxY6r1WoaNiw72OedxpO9rLOR01bM+OMMB6/eYO5TLS2yGdJUJv6vETvPJ/LGmsP8+nJnXB3t57FXnhVchcuaslHSlthzznv//CjFht8lRqsTQqiFEIeBeOBPKeV/pZR7QQgRLoQIT0ioWKIeBQUF09h2Ko6v/73I4A41+F9zy2dYLQsvV0dm9G/Oufg03vjhCDoLhdW3JoXjhNljmBebOhYhxFYhxPESfkyOqCilzJNStkCfh6WdEKJJKeUqlOhLQUGhfERez+CNNUdoVN2T8X2sk57ZGF3r+/H+Qw3ZfCKWyRtO3hHOJR97C5kPNh4Kk1L2KO2cECIuP3e9EKI6+h5JWXXdEEJsB3qjj1mmoKBgZa6l5zD823B0Osk8G82rlMawLrWJvpHF0l2XiE/JZvoTzXA3Q5IxW2MPGyKLYs9DYeu4mdflWfTJxW5BCOEnhPAyfHZBv6HktNU0VFBQKOCv03E8MOtfLiSkMe/pVjadVykJIQQfPNyQ9x9qyMbjMfSY+Q/rjkRbLOPk3Yw9u+tpwBohxDDgKtAfQAgRCCyWUj4EVAe+FUKo0TvJNVLKstP6KSgomI3Y5Cx+OhjJjnMJ7L14jQYBHiwf2o5GgZZPW1ARhBCMuDeUNrW8GffzMV5dfYjPt5zhkZZBdG/gT8PqHmjU9vy+fXtgt45FSpkEdC/heDT6jJFIKY8Clg08pKCgUCpJ6dl89scZGgR48EbP+rzYNRQnjf0Mf5VGyxrebHj1HjYfj+XbPZeZve0cX2w9xwON/Fk4xDIbcO8m7HaDpCWpWrWqrFWrlq3VsAsuX76MYgs9ii1ucvDgQVxcXGythl3g4uKitAsDBw4cQEppdHOZ3fZYLEmtWrWsttvc3rHmznt7R7HFTRRb3ESxxU2EEBnGS9n35L2CgoKCwm2I4lgUFBQUKkHk9QyW7LxE+OVrtlbFbrgrh8IUFBQUzMHVpAz6ztvJjQx9jhV7yZtia+yixyKE6C2EOCOEOG8IOFn0fAMhxB4hRLYQ4q0i57yEEGuFEKeFEKeEEB2tp7mCgsLdzCcbT5GXJ1k3ujN9mlVn6qbTBdkj72Zs7lgMe1DmAQ8CjYCnhBBF40BcA14FZpRQxWxgs5SyAdAcOGVBdRUUFBQAiLiWweYTsTzbqRbNgr2Y2b85DQI8eP+X46Rl59paPZtic8cCtAPOSykvSilzgO/RRzYuQEoZL6XcD2gLHxdCeAL3AksM5XKklDeso7aCgsLdzObjsQAMMMTocnZQM+XRpsSmZDFn2zlbqmZz7MGxBAGFw3JGGo6ZQiiQACwTQhwSQiwWQpQYR0KJbqygoGBOtpyMpVF1T2r4uhYca13Tm8daBvPN7svEJGfaUDvbYg+OpaTNNqbu2tQArYD5UsqWQDr6hGDFK1SiGysoKJiJLG0eh67e4N76xZ8lr/Woh5SSOdvO20Az+8AeHEskUDjeczAQXY5rIwvlYFmL3tEoKCgoWIyjkcnk6iRtanoXOxfi48qgdjVYEx5B5HWT9hPecdiDY9kP1BNC1BZCOAID0Uc2NoqUMhaIEELkJ4vuDpy0jJoKCgoKesKv6PestCrBsQC8YFhy/O3uy9ZSya6wuWORUuYCo4E/0K/oWiOlPCGEGCmEGAkghAgQQkQCbwDjhRCRhol7gFeAlUKIo0AL4BPr34WCgsLdxKGrNwit6lZqKuAgLxf6NK3O9/siSM3SlljmTsYuNkhKKTcCG4scW1Docyz6IbKSrj0MKOFIFRQUrMapmBRa1ii5t5LPiHtCWXckmrUHInm+c20raWYf2LzHoqCgoHA7kZadS+T1TML83css1zS4Cs2Dq/DD/gjutijyimNRUFBQKAdn41IBCAswnsysf5sQTsemciI6xdJq2RWKY1FQUFAoB2djDY7F38No2f81D8RJo+KH/RFGy95JKI5FQUFBoRycjk3F1VFNsLfxRGhVXBx4sEkAvx2OIkubZwXt7APFsSgoKCiUg7NxqdTz90ClMppIEdCHfEnJyuWPE7EW1sx+sItVYXcCqampZGRk4O/vX+m64uLiyM3NxdPTEw8P493tO4UrV67w06a/OX/2NEIIpJTodDp8/arRpn0nHuneESFM+zLbM79s28PxwwdwdffA29uHgQ93x9XV1fiFCnbBmdhUujesZnL5DqG+hPi48MP+CPq1MDVa1e2N4lgqyZfLVnPgv934+lXDycmZuJgo/KsH8fboEXh5eSGlLHhIAsUejDqdjpMnT/LTpm3EREUCUNXPHwdHB1KTk0lLS0UIQXCNWgx76jECAgKsfo+VRUpJRkYGGo2GnJwc0tPTyc3VR3/dfuAkB/buIjU1heCQmjRr1ZYHHn4ElUpVcG1SYgL7d+/ghd/W4u7hyVsvjyAo6Pb6gsbHx/P5/KUkxMdRN6wh7TrfS2ZGBkmJ8YydNJXMjAz8qwfy9KMP4e3tTVZWFunp6SQnJ3P4fBRZmZnk5eXi5x9Azw7NqV69eoGN7B2dTkdycjKJiYkkJiay/8QFYqIjSUm+S4+UNQAAH/xJREFUGS+28Pei6Aqqwt+f/PNqtZrA4BrUC2vIoz27oFarLX8jQFJaNknpOdQ3YX4lH5VKMKB1CDP/PMvVpIxbYovdqZjNsQghVMBYKeVdsUExOzub18dOoFGzFrwzccot56IjI/h45lwyMzLKXmYoJQhB3fph3HN/TwICg0p9I7984TxzFi0nKTEetVpD3bAG9OtxD6GhoXb5Fh8dHc2ilWuJjoxApVLh4uqKLi8PB0dHXF3dUKnVCCGo5l+dQc+/gJt7yV9UIQRV/arxYL/HebDf49y4do05i5aTmBBH53u78Vz/h236gI2KimLt5u14eHjSvX0zQkJCCvS5du0a85f/wMVzZ/H28eXBfo/jXz3wluvrhjWkfeeuAMRGR/HrnzvISE/D2cUFJ2cXPDw88fL2wTnQBbVaTXxcLAu/W0t8bHRB23L38CS0Xn1a1gvh2KVY2jasRc2aNfH29rZZ28jMzGTC1Jmk3LiBg6MDHp5V8PbxxdvHlxq1Q2nf5V48q3hVWL/c3FxioyM5dfwYr7z9HlJKPDyr0L7TPTzSo3MxR6PValn52x/s272Dth278GSfivUSz8WnAVCvHI4F4Ik2wczaepY14RG81SvM+AW3OcKc66uFEH9JKbtV4Lre6POqqIHFUsppRc4Lw/mHgAzgOSnlQcO514Hh6ANXHgOel1JmlSWvTZs2Mjw8vLxqFhAREcHYCZN5/qVXCAqpWeF6Kkpubi7nz5zi5LEjXL18EZVKRZf7ejCob89yf1HbtGlDZWyh1WpZtOoXThw9VPCgy9Vq8QuoTpeu3akZaplserm5uezduZ29O/8hOKQmH777WqUdTHlssXX/Cb5ZMJeaoXVp1qotGWmpXL18ibiYKHQ6HRKJu7snXXv0IrSeZR8kqSnJXDx3lpTkG7i5u3M9KZGYqEiSb1wnJyeH2nXrMWLQ4+Uapq1Mu1j3917WrlrO8yNfoXpQifuaLcKN69fYt2sHZ04dLxgpEEKQm5uLRqOhZZv2tGzXgaMHw9m87hdWLltoUk+nsC2+23uF8b8eZ9fYbgR5GZ+8L8xzy/ZxOiaVne/ej0Z9e/Q2iyKEyJBSlhhB/pZyZnYsM4EUYLKUUmfiNWrgLNATfVDJ/cBTUsqThco8hD50y0NAe2C2lLK9ECII2Ak0klJmCiHWABullN+UJbMijiUzM5Pv129l1/Zt+PkH8NjAwbh7GF/Hbg20Wi3/btvCkQP7WDj383I5l4o+QDIyMli8+hd2bd/Gw48NoGXbDjbrOZw8dpgNv/zIrGmT8fLyqnA9ptpCp9MxaOiLvP3Bxzg5O1dYnrW4eP4sa5YvYcgLL/NA+6YmXVPRdpGWlsao195m7EfT7Hqo7uypExzYt4dJ775utGxhW3y47gRrwiM4MalXuV/itp6MY/jycGb0b84Tra3ncM2JqY7F3P/5EPRBJKOFEL8JISYLIfobucZooi/D38ulnr2AlxCiuuGcBnARQmgAV0yPjGwSsbGxjBj9Ou9OnIIuL48xYycwZMQou3EqAA4ODnTv3YfuvR9m/JSSkmyaFyklo98ah3/1QD6Y+jmt23ey6UOkUdMWPP/SGN77aBoz5i+1uLxTp07Rul3H28KpAITWrc/r703imwVzLb4D/L2PpjH8lTfs2qkA1G/YmLjo8j8qLiSkUcfPvUJDeN0bVqNxoCdz/zpHdu6dvfTYrP99KeUAKWVDoCYwCTiP3nGUhSmJvkosI6WMQp+u+CoQAyRLKbeUJKQiib627T/Bm+Mm8NLr7/LS6+/QqWs3NBr7Xe/QvHVbkhLjLSojNzeXYaPGcG+3B2jaorXdPECq+lVj9Fvvoc3RMnbSNDIyLBeufMP2vdRv2Nhi9VsCJ2dnAgKDSU9Pt5iMS5cu4ermRjX/22OBiVCJcjvac3Fp1KtWdiiXUuUJwTu9G3AlKYMv/7qzc7VY5KkgpcyWUh6UUn4rpXzbSHFTEn2VWEYI4Y2+N1MbCATchBCDS9Gp3Im+vl++hHcmTsHlNloKqhKWfdAfP36cpi1a067TPRaVU1H69X+Krt17MfSlV/n7wGmLyLhy8Ty16tSzSN2W5Ma1JNzdK/ZQNIXPv1rEo08+bbH6LUF5eh6pWVpiU7KoU0HHAtC1vh+Ptwrmq+0X+Pu0ZV8CbYk9vG6akuirtDI9gEtSygQppRb4GehkLsUcNA44OjmZqzqLk56WSm5erkVlZGZm4lrKCi57oWZoHd6e8DHzPp9qsaGf2y2oYPjeXdRv1MTiclxcbp+XsPIOZ11I0Pf26lbCsQBM6teYhtU9GLXyIOGXr1WqLnvFHhyLKYm+1gFDhJ4O6Ie8YtAPgXUQQrgaVo51R5/TxSx4+fgSGx1lruoszqpli5jwjvHJyMrg4+PDhTNmM7HFcHFxxdvH1yLLbVu07cCu7VvNXq+lyM7KYtvm9bw7erhF5TRq2oKTx49YVIa5SE1JLnWJe2mcNyw1rqxjcXfSsOy5dgRUceb5Zfs5GnnD+EW3GTZ3LKYk+kKfq+Ui+jmbRcAow7X/oU9HfBD9UmMVsNBcur05ahi///SDuaqzKCePHsbN3YPgYMuuNgkLC8PByYnzZy0zzJRPRkY6Z0+dIHzvLg7u38uVixfIyyvfhKejo2V6m8Of7MuBvbuJibo9AguuWPwVzwx/yeJ7Wi5fPE/N2pZZXm5u/t6yiS73dS/XNefj03BQC2r6VL5X5ufhxMrh7ani6sAzS/Zx8g6LfmwXM9EmJPqSwMulXDsRmGgJvXx9fUlNSS5YE19ePl65lQv/bSVXm4MwxBWSOonGwZHqDVry8v86oNXmkppyg2/+PEhmynW0WRlIqcPZ3QuPqgG80b8bfkYmQ3Oys/llzSqWL/qqQvdZXj58ZwwjRr9RbGNoWUxZtZXYc0fJzcnC2b0KLp4+DLmvCa5ubqjVar7efJCkiHOkJenjKWkcnfH0D8bZ3QupyyMjeS/JcZHocrX4htTl07deMDpMqdOZtOK93AghmP/Fpzz74mjGfTS9Uru+8/LyiLp6hc9Xb+RGzBX9ptlbphQlGidnqtZswPjn+5Z7NeKZk8fwrOJl8jLjypCellruxS1pqSmM//I7rkVeLNgwrJ9iLfpbj0qlwdHVHQx7VJzcPPHwDeC9Z3qXqwcSceUSb740tFy6no9PpXZVN7PtQQn0cmH1iA4M+HoPzyz5j+9f6FDujZf2il04FnsmMDiEhLhYqgVUN14Y/dj7e/NWc+XwLryq16RZ76dwdLl12bc2O5PYs0f4ZPFa1BoHnFzd8fANoFrtRji4uCKEiqy0ZFITovlg9lLSkuJw8fRm3sfjSnyYpqWl4uvnZ7UVWg4ODrRo046D+/fSqm2HMssmxMfx2vsf4R0USlDDVjg4u5KVmkxGchKLNu5Bm5WBLi8PVy9fajTvhLuPv1EnnnD5NM+OGkODLn2Y+Pz/SiwTFxONmwUnqp2cnLi3W0+OHtxPyzJsEBMVwYRZi9BmZd48WPhZKcDTLwi/2g2p074nKlVxJ5WTmU7ilTO89uFn5GSmoXF0pmbzzkx8/mGjev66ZhXfLJhbzrurGB+8/RpTZs5l1BvvGi07be1Ojm9di4OTMzVb3EOddt1LvPei6PJyyc5IAymRUpKdnkJKQjSvfDCdnIxUnNw8mPzGiwQGh5RZT15eXrlfCKY/3oxr6TnlusYYIT6urBzenicX7mXwkv/4ZVRnAsu58dIeURyLERo2bcHxI4foZsSxZGSk8+YnX5IcF0FQw9Z0emoMopQHvYOTCyFNy34gO7q44ekXSFAjfdbllPgonhv9JqsWfVmsrLu7BzeuWXcS8LFeXfn2x3VlOpZPf97NofXf0n7Ayzg63xw+cK3ii09waIVl+9VqQNWaYexcMQPds32KOdS8vDwWz5vFoi9nVViGKQwd0Jc3x39UqmN5ecJnaLPSadC1H87uVSosx9HFjcAGrQhs0AoAbVYmlw/+w6BhL9Gk55O8N/C+Eq+7euki9cIaWi2Olr+/Py6ursRGRxEQWHost/Ff/8SVwztp028YDs7le4iq1BpcPG5ugnWt4oN3YC1qNtev2clKS2binKVkp6ey5IupODg4FKtDp9ORk1N+B+Hr7oSvu/mHV0P93FkxrB395+/h+WX7+fGljng6F9f7dsLmcyz2zhO97mX3v39zLSmx1DJvz1zKC6+PI7hxOzoPeo1aLe8p1alUFM9qQWgcHIsdl1Ly5YxPeGbES2aVZ4zAwEBOnzhGelpqqWWO/vE9HZ4cfYtTMRdCCGo068SEr9cWO/fdkgU8MehZi0cMdnZ2Jjsrq8QVYjnZ2aRfi6f5g4Mr5VRKwsHZhXqdetNh4Kuc2LqGG9dLfqn4+8+NvDDY2P5k8/L6yKFs+LX4/yQfKSXn926hff9R5XYqpuDsXoUWDz5NWOeHGPnupBLLbN20nvt7Pmh22ZWhQYAnC55pzYWENEZ9dxBtnmWGca2F4liMoNFomP/FZyz9anaJ56MjI0iKvEDHga/iFVB297sySCnJzckudnzF4vl07dGbHm0tv5S0MBqNhtmfTmH29Mmlvv2pHRxwcLJct75KQAipSbfmuDhyYD8uLi482sNsq87LpFpAdRLj44odV2s0mBjVqMKoNQ407jGAD2YvK3ZOSsm1xESrR8MOCAgg+fp14uNKzj2SnpaKm3c1iy8k8Kpeg5zM4ptkpZQc3LeHp/s9YFH5FaFz3apMe7wZO88nMu7nY7fdkvbCKI7FBDw9PWnSohW7//272LnZP22nelgLi+twZsd6arW695Zjv/ywksDgGgx+pJfF5ZdE1apVGTJiFHM+nWyxifKySEuKxdWrasHfOdnZbPxtLRPfGWM1Hbr1eoi1K78pdlytViOEipzMNIvKV6k15OVqix2Pj43Br1rlcwNVhC+mT+br2Z+V6Fzc3D3ISku2uA7ZGamoNcWHk/76YyP39extcfkV5YnWwYzpXo+1ByIZ/+tx8nS3p3NRHIuJvDZiCHt3bOfKpQu3HH//uf8RdWK/RWXHnj9GVnoKn4x6suDYts0bAHhl6CCLyjbGAx2a0evhR1ixeH6xc7q8PKQFHU7UyXA+GHozrNzyxV8xeNhIq4aK79G2MT5V/Th76kSxcw269uXkX79YVP653RuZ+PKQW47lZGfzzddzeXXEkFKusiyurq4s+vILFs6ZQXxszC3nhBA4OLmQk2W5kDsAB9Z9w4wJb95yLCc7m/92/ctzTxhf9GBLXutRjxe7hrLyv6uMWB5u9gUD1kBxLCYihOCrWZ+yYvECriclFRz3rOKFSuNASkJMGVdXnBsxV7gU/jcLp08oOBa+dxeRVy/z3mvWnVcpjScfup8qXt78t+vfW47XaNqBSwf/sYjMhMtncK1SFVdX/Yq7MyeP4+LqRq+OzS0irywmvPUKf/2xga0bb93XO7a/PuzN9ejLFpF7KXw7Hn6Bt+R4ibx6mc8+Gs+QEaMIDAws42rLUuBc5s4kLubWQBr1Oz/IiW0/WUz2sT9/JKRJO6r63ZrlcdWyRQx6foTF5JoLIQTjHmzI5H6N2XkukV5f/MumYzG31dCY4ljKgYODA/O/+Ix5M6cWZEAEmD/1Aw6t/9bsb2EpCTEc2/oj38y7GQr/5LEj7NnxD1MnjDWrrMoybsyL7Pj7T6ILbRr8+KUniT51kKx0827+Sk2K5cyO9cz7WG+D3Nxc1q76lo/GWjbqQGk4ODjw5WdT2LLht2IT6V9NHc/xP38gK9W8u6sv7NtK+vV4vnhvdMGx7X9u4rc1q/h24Zf0bGfdObeScHFxKXAu6ek3hwTfG3gfUqczu8OVUnJk0yrcfaox7ZVbY5ZFRVxFq82hd6eWZpVpSZ7pWIvfRnfG182Rl1YeZODCvRy6et3WapmEXTgWIURvIcQZIcR5IUSxJ6YhlMscw/mjQohWpl5rbjw9PXl66IvM+uTDAufi5OzM/M+nsmf1HNKumSewXOLVcxzdvIrl82cXbDrbs2M72zZv4MsZn9hd1kghBF99Pp0fv/uG3f/+VXBswedT+W/NPLM53dhzRzm6eTXL5s1CpVKRnp7Gp5PeY/DQkVZbVlsa3y1dyPKF8/j9px8K3i4dHBxYNGcG4b8s5FrUxUrL0GZlsv+nBahUGuZ/8h6g32Q4e9pH5Gq1zJ81HUfH4qsHbYWLiwuPDBjEru1/3XJ8wfQJnPjrZ5IiL5RyZfnIycpgz/dz8KvdkBlvDSt2/rslC/h4/DtmkWVNGlb3ZP0rXfj4kSaci0/j0a92M2DBHraciCXXjleOmTXRV4UUqFyiL6PXlkRlM0gC/P7PPtasWMbI197Gx1c/gZyZkcHIdybiWsWHhl37oS5hebAxdHm5nPpnHVlpySya8RFqtZqcnByWzp9NYHAN3nl5uFmdSmUzSJbElFnzibhyiQf69KNB46YkJSYw+p3x1GrVlZAmxrIolExG8jWO/bkGd59qfDnp7QIbaLVaQjwEPj4+ldbbXLZY/stmtqz/laeee6Egg2Zubi6jxn2MTpdHkx79cSjnEmyp03H16G4ij//HF1MmFGzY3bfrX7Zv3cyMjydSrVo1I7WYjjnbhVarZeyHn6DVanlyyFB8q+qji+fl5THstffwCQ6lTrvyhVfJR0rJlcM7iTj+H3OnTSpxwcLfWzYhVIKRFVx6bYnvSEVIy87lh/0RLN15iagbmfh7OvFE62AGtAmhpq/R3FtmwSYZJCuCEKIj8KGUspfh73EAUsqphcp8DWyXUq42/H0GuA+oZezakjCHYwG4fv06H0z5DG8fX/oPfq5gM9bk7/7kzM4NaBwcCWnWEf/Qxkb3taTER3Fh/zYyk69Tr1NvPny+DwDHDh/gtx+/Z/CwkTzY2fzdeEt9aTIzM5m7dBX33N+DKl7eSCl5e+Ziok8fIqBeU2q3vq/EVTuF0enyiD13lKtHduPg5MKsD9/Bs8qtGSIbBZov4Zo5bZGTk8P7Uz5Dm6Nl4LPDCkKxxEZHMvaT2ajUamq37YZPUNkbRVPio7h6ZBepiTGENO3AJ6MHIYQgMzODZV/NpkbtOrw72rwvG2CZdpGUlMSUmXPR5mp5tlCyvLdmLiXyxH4C6jahRvNOOLkaD2uSm5PNhf1/kXDpFCFN2zN9zJACGyTfuM6m334mPi4GIQRZWVksmfdFhfW2F8eST26ejq2n4lkTHsH2M/HoJLSr7cP9YdVoVcOLhoGeFttgeTs5lieA3lLK4Ya/nwHaSylHFyqzHpgmpdxp+Hsb8C56x1LmtYXqeAF4AaBGjRqtr1y5YrZ7WP/vftasWEa//k/RqOnNyePMjAzGz1tFwuVTgD62kZQSByfnW/alSJ0Od59qTHljeEHv59jhA2zdtJ7gGjUZ/8Zoi4VrsdaXpnCQvQmLfubKkd1InQ6VWo3G0QkHZ1eklOSkp5KXZ1g+KyX+dZow5ZVniu2gNqdDyccStoiKiuKjz77Ax9ePJwY9W5B1MjUlmfGzv+F69GWESoVQqVBrHPU52nOyAX3IEnefanzw0tMEBN4MLrpr+zZ2bt/Ksy+MpkfbRmbVNx9Ltov4+HjGT56Oj29VBg8fiUajIScnh09W/knEsT1kp6ei0mhwreKLSxVfNA6O6PLy0GZnkn4tjpzMdFRqDbVa3cukofqQPokJ8fyx/lcSYmNxc3fn5WHPEBoaWjAkWRnHa2+OpTAxyZn8dCCS34/EcCbu5mbl6lWcCfF2JdjbhSBvF3zcHHFz0uBu+HFzUuOgVqFWCTSq/N8CjVoQ7F16b/p2ciz9gV5FnEM7KeUrhcpsAKYWcSzvAKHGri0Jc/VYCqPT6Zg4bRaJCXGo1Rq69+5DWKMmxRq0lJKszEx9AD1n51Ib/P49OxnwUDecLZz+1tpfmqJRXKWUZGVlkp6Whkqlwt3Ds9Q5Aks4k8JY0hYbdx7i5+9X4FfNn34DBhULJimlJCc7G53U4ezsUmK7OHPyOOvWrqZVu468OnSQRefZrNEuNu08yHdLFxBaL4w+jzyBl/fN4UwpJQnxcVxLTCAnJxu1WoOrmxtV/arh4amPZJCVmclfWzZy8tgRfQbRYc9YJLq3PTuWwiSkZnM8OpmT0SlciE8j8kYmUdcziUnOxNTtME4aFWc+Lj0qwe3kWKw+FCaESABK6rJUBUqP3WI+rCXHFFmt0KcdqGw95sDW9i/LFrbWzdqyTGkX9qKrpeW0BvJXn2gAy2bTs285LlJK48Mn0hAl1FY/hhu7iD69sCNwBGhcpEwfYBP6uLAdgH2mXltOXcKtdM9WkWNOWdbQ2Z7tb8+62UrW7aSrueTYix72Lsfm0Y2llLlCiPxEX2pgqTQk+jKcX4A+V8tD6BN9ZQDPl3WtDW5DQUFBQcGAzR0LVDrRV7FrFRQUFBRsh11skLQjzJbW2E7kmFOWNXS2Z/vbs262knU76WouOfaih13LsfnkvYKCgoLCnYXSY1FQUFBQMCuKY1FQUFBQMCt3tWMRQvgIIf4UQpwz/PYuoYyzEGKfEOKIEOKEEKLkfKcl11/h4JoVuBdjsvoZZBwWQoQLIboYqc9itrGGXSpjD0u3CxP1M0vbMGe7uJPtUsLfxexSmlxjOpvr/oQQrxtselwIsVoIUeruaRPkNBBC7BFCZAsh3ipyzksIsVYIcVoIcUro9xqWD2uslbbXH+BTYKzh81hgegllBOBu+OwA/Ad0MKFuNXABfXSA/D02jYqUeYhb9+f8V8H7MEWWOzfn1JoBp21hG2vYpbL2sGS7sGbbMHe7uMPtcrXI322K2CWiJLmm6GyO+wOCgEvoNygCrAGeq4ScakBbYArwVpFz3wLDDZ8dAa/y2viu7rEA/dAbEcPvR4oWkHryk0k4GH5MWfHQDjgvpbwopcwBvjfIKyp/uUHGXsBLCFG9AvdhVJaUMk0aWgrgZsI9WMo21rBLZe1hyXZhkn6Yp22Yu13cqXbZDWQVkduziF3cS5Fris7muj8N4CKE0ACuQDQlY8r/PV5KuR+4Ja+1EMITuBdYYiiXI6UsdzKhu92x+EspYwAMv0uMOy6EUAshDgPxwJ9Syv9MqDsI/VtOPpGGY+UtYwom1SOEeFQIcRrYAAw1UqelbGMNu1TWHpZsF6bqZ462Ye52cafaRcutD9hIIKiIXU6UIrc8+lT4/qSUUcAM9D2rGCBZSrmlEnJKIxRIAJYJIQ4JIRYLIcodk/+OdyxCiK2GMcmiP6W9VRRDSpknpWwBBAPthBCmpOcrKUJg0Tc3U8qYgkn1SCl/kVI2QP+mOdlGtrGGXYxeL4TYCkxGHyfpBvCnldqFSfqZWMYccoq2i/U2+r6Yqq8l7FJinUXs0qAUueXRp8L3Z5jP6oc+fFUg4CaEGFwJOaWhQR8nbr6UsiWQjn7Ys1zYxc57SyKl7FHaOSFEnBCiupQyxtDdLDP9o5TyhhBiO9AbOG5EdCQQUujvYIp3XU0pYwrlqkdK+a8Qog7QVkpZYmA/C9rGGnYxen3RdiGEuATcJ6VMtHC7MEk/E8uYQ04BhnaRgsEORc/fwXbRoB+yK7FOg10cgUbAziJlHMuhT2XurwdwSUqZACCE+BnoBHxXQTmlEQlEFuplrqUCjuWO77EYYR3wrOHzs8BvRQsIIfyEEF6Gzy7o/8GnTah7P1BPCFHb0CgHGuQVlT/EsBKkA/rubUwF7sOoLCFEXSH0cdYNK00cgaQy6rSUbaxhl8raw5LtwiT9ME/bMHe7uFPt0gn93EVhuYeK2CUHeLQEuabobI77uwp0EEK4GvTqDpyqhJwSkVLGAhFCiDDDoe5AmRl5S6vorv0BfIFtwDnDbx/D8UBgo7y5IuQQcBT9W9eEctT/EPrUyReA9w3HRgIjDZ8FMM9w/hjQphL3YkzWu+jHiQ8De4AutrKNNexSGXtYul1Ys22Ys13cyXYx/B2Hfn7hfYNdotHPVewBupQmtySdLXF/wCT0Tvo4sAJwqoScAPS9kxT0Q8GRgKfhXAsg3PA//BXwLq99lZAuCgoKCgpm5W4fClNQUFBQMDOKY1FQUFBQMCuKY1FQUFBQMCuKY1FQUFBQMCuKY1FQUFBQMCuKY1FQUFBQMCuKY7EiQh9DabbQh74+JoQItbVOtkKxxU0UW9xEscWdgeJYrMs44KKUsjEwBxhlY31siWKLmyi2uIliizIQQqhtrYMp3PGxwuwFoY8Q+qiUsrXh0CWgjw1VshmKLW6i2OImii1KRgjxI/oIAC3RRzz42LYaGUdxLNajBxAi9OHEAXyArTbUx5YotriJYoubKLYomabAKSnl/bZWxFSUoTDr0QJ93KQWUh9SfAtwWAjhJoT4VgixSAjxtI11tBal2SJUCLFECLHWxvpZk9Js8YihTfwmhHjAxjpai9Js0VAIsUDo0+W+ZGMdrYrQpx/2AT6ytS7lQXEs1sMbyAAQ+gxwDwC/A48Ba6WUI4C+tlPPqpRoC6nPeDfMpppZn9Js8auhTTwHPGk79axKabY4JaUcCQxAnzL4bqIx+vTEubZWpDwojsV6nEWfwxrgdWCDlPIS+lwJ+dne8myhmA0ozRZ3I8ZsMR59tNu7gVJtIYToiz4XyjYb6WYrmqKPMnxboTgW67EaaCWEOI8+tPgbhuOR6J0L3D3/j9JscTdSoi0M+TimA5uklAdtqaAVKbVdSCnXSSk7AXfLcHE+t6VjUcLm2xjDSpgvgSxgp5RypY1VshlCCF9gCtATWCylnGpjlWyGEOJV9Mm09gOHpZQLbKySzRBC3Id+yNgJOCqlvFt6cLctimNRUFBQUDArd8vQi4KCgoKClVAci4KCgoKCWVEci4KCgoKCWVEci4KCgoKCWVEci4KCgoKCWVEci4KCgoKCWVEci4KCgoKCWVEci4KCgoKCWVH9cSJW/nEiVtklqaCgoKBgFpQei4KCgoKCWfk/uOfAZMyYMPsAAAAASUVORK5CYII=\n", - "text/plain": [ - "
" - ] - }, - "metadata": { - "needs_background": "light" - }, - "output_type": "display_data" - } - ], + "execution_count": null, + "id": "1ec02909", + "metadata": {}, + "outputs": [], "source": [ "try:\n", " from anesthetic import NestedSamples\n", @@ -213,34 +212,9 @@ "\n", " print(\"Install anesthetic or getdist for for plotting examples\")" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.8.2" - } - }, + "metadata": {}, "nbformat": 4, - "nbformat_minor": 2 + "nbformat_minor": 5 } From 8954de3c69d719763b2344553ef6adb4a24f002a Mon Sep 17 00:00:00 2001 From: Ormorod Date: Thu, 22 Sep 2022 12:22:20 +0100 Subject: [PATCH 32/39] correct docstring to clarify cluster function must number clusters from 0 --- pypolychord/polychord.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pypolychord/polychord.py b/pypolychord/polychord.py index eea20fba..20cc19bc 100644 --- a/pypolychord/polychord.py +++ b/pypolychord/polychord.py @@ -118,7 +118,7 @@ def run_polychord(loglikelihood, nDims, nDerived, settings, Returns ------- cluster_list: array-like - cluster labels. Must start from 1. All clusters must have at least + cluster labels. Must start from 0. All clusters must have at least one point, so that max(cluster_list) gives the number of clusters found. Length nPoints. From ebf3b91d8310cc3763d7d6db8635c16279206f22 Mon Sep 17 00:00:00 2001 From: Ormorod Date: Thu, 22 Sep 2022 12:24:27 +0100 Subject: [PATCH 33/39] further correction to clustering in docstring --- pypolychord/polychord.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pypolychord/polychord.py b/pypolychord/polychord.py index 20cc19bc..47db9cea 100644 --- a/pypolychord/polychord.py +++ b/pypolychord/polychord.py @@ -119,7 +119,7 @@ def run_polychord(loglikelihood, nDims, nDerived, settings, ------- cluster_list: array-like cluster labels. Must start from 0. All clusters must have at least - one point, so that max(cluster_list) gives the number of clusters + one point, so that max(cluster_list)-1 gives the number of clusters found. Length nPoints. From e1102fbd903a2b999e653324fe3b35e0376b0072 Mon Sep 17 00:00:00 2001 From: Ormorod Date: Thu, 22 Sep 2022 12:34:46 +0100 Subject: [PATCH 34/39] change cluster function example comments to docstring --- run_pypolychord.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/run_pypolychord.py b/run_pypolychord.py index e03c818b..49115948 100755 --- a/run_pypolychord.py +++ b/run_pypolychord.py @@ -42,16 +42,17 @@ def dumper(live, dead, logweights, logZ, logZerr): #| Optional cluster function allow user-defined clustering def cluster(points): + """ + + - clusters should be an array of cluster labels for each point + - each cluster should have at least one point + - thus max(clusters)+1 should be the number of clusters + - i.e. clusters are 0-indexed + - work with the above numpy integer array + """ npoints = points.shape[0] clusters = np.full(npoints, -1, dtype=int) - # - # - clusters should be an array of cluster labels for each point - # - each cluster should have at least one point - # - thus max(clusters)+1 should be the number of clusters - # - i.e. clusters are 0-indexed - # - work with the above numpy integer array - return clusters #| Initialise the settings From 4d3adcc164a68076d0109dc9480a1e08239e95d5 Mon Sep 17 00:00:00 2001 From: Ormorod Date: Thu, 22 Sep 2022 12:45:48 +0100 Subject: [PATCH 35/39] remove unnecessary np.sqrt import --- run_pypolychord.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/run_pypolychord.py b/run_pypolychord.py index 49115948..c33d54ed 100755 --- a/run_pypolychord.py +++ b/run_pypolychord.py @@ -1,4 +1,4 @@ -from numpy import pi, log, sqrt +from numpy import pi, log import numpy as np import pypolychord from pypolychord.settings import PolyChordSettings From 70628a5bb947f23bb248dde3ebeeaf650efa0b61 Mon Sep 17 00:00:00 2001 From: Adam Ormondroyd <52655393+Ormorod@users.noreply.github.com> Date: Tue, 7 Feb 2023 18:49:02 +0000 Subject: [PATCH 36/39] Correct root=root -> root=0 in make_resume_file (#92) * neighbor spelling (assuming we're sticking to American spelling!) * change NameError catch to check whether mpy4py was installed * set root=0 rather than erroneous root=root * Revert "neighbor spelling (assuming we're sticking to American spelling!)" This already appears in a different PR This reverts commit f1f52edbd5b8a9593342864dfd1d0f25ac0c8c49. --- pypolychord/polychord.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/pypolychord/polychord.py b/pypolychord/polychord.py index 580fba94..1d497d58 100644 --- a/pypolychord/polychord.py +++ b/pypolychord/polychord.py @@ -232,6 +232,7 @@ def make_resume_file(settings, loglikelihood, prior): rank = comm.Get_rank() size = comm.Get_size() except ImportError: + MPI = False rank = 0 size = 1 @@ -245,17 +246,17 @@ def make_resume_file(settings, loglikelihood, prior): nDerived = len(derived) lives.append(np.concatenate([cube,theta,derived,[logL_birth, logL]])) - try: + if MPI: sendbuf = np.array(lives).flatten() sendcounts = np.array(comm.gather(len(sendbuf))) if rank == 0: recvbuf = np.empty(sum(sendcounts), dtype=int) else: recvbuf = None - comm.Gatherv(sendbuf=sendbuf, recvbuf=(recvbuf, sendcounts), root=root) + comm.Gatherv(sendbuf=sendbuf, recvbuf=(recvbuf, sendcounts), root=0) lives = np.reshape(sendbuf, (len(settings.cube_samples), len(lives[0]))) - except NameError: + else: lives = np.array(lives) if rank == 0: From 6cf0fc4a9bc462d58673b4687d2e7cb879ba4a59 Mon Sep 17 00:00:00 2001 From: Will Handley Date: Thu, 16 Feb 2023 18:33:28 +0000 Subject: [PATCH 37/39] PWD->CURDIR --- Makefile | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Makefile b/Makefile index 8e871ce0..adca6c15 100644 --- a/Makefile +++ b/Makefile @@ -5,13 +5,13 @@ EXAMPLES = gaussian pyramidal rastrigin twin_gaussian random_gaussian himmelblau PROGRAMS = polychord_fortran polychord_CC polychord_CC_ini # Directories -SRC_DIR = $(PWD)/src +SRC_DIR = $(CURDIR)/src DRIVERS_DIR = $(SRC_DIR)/drivers POLYCHORD_DIR = $(SRC_DIR)/polychord -LIKELIHOOD_DIR = $(PWD)/likelihoods +LIKELIHOOD_DIR = $(CURDIR)/likelihoods EXAMPLES_DIR = $(LIKELIHOOD_DIR)/examples -BIN_DIR = $(PWD)/bin -LIB_DIR = $(PWD)/lib +BIN_DIR = $(CURDIR)/bin +LIB_DIR = $(CURDIR)/lib export DRIVERS_DIR POLYCHORD_DIR PYPOLYCHORD_DIR LIKELIHOOD_DIR EXAMPLES_DIR BIN_DIR LIB_DIR # Whether to use MPI From 3084a819f3f37aead107e59967013cbbfd72776e Mon Sep 17 00:00:00 2001 From: Will Handley Date: Thu, 16 Feb 2023 21:32:08 +0000 Subject: [PATCH 38/39] Updated PWD to CURDIR in setup.py --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 17b939eb..23ff9384 100644 --- a/setup.py +++ b/setup.py @@ -92,7 +92,7 @@ def run(self): env["DEBUG"] = "1" BASE_PATH = os.path.dirname(os.path.abspath(__file__)) - env["PWD"] = BASE_PATH + env["CURDIR"] = BASE_PATH env.update({k : os.environ[k] for k in ["CC", "CXX", "FC"] if k in os.environ}) subprocess.check_call(["make", "-e", "libchord.so"], env=env, cwd=BASE_PATH) if not os.path.isdir("pypolychord/lib/"): From 185437268b0ef64cfc6b72fa4b023bee10d7c988 Mon Sep 17 00:00:00 2001 From: Adam Ormondroyd <52655393+Ormorod@users.noreply.github.com> Date: Mon, 20 Feb 2023 12:12:13 +0000 Subject: [PATCH 39/39] Update python3.6 CI (#96) * neighbor spelling (assuming we're sticking to American spelling!) * change to checkoutv4 * change to checkout@v3 * change to pythonv4 * change to ubuntu-20.04 * return to latest ubuntu for 3.7 and 3.8, and separately include 3.6 case * include pythons up to 3.11 * enclose versions in quotes * revert spelling changes --- .github/workflows/CI.yml | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 4284765e..66f875f2 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -16,12 +16,15 @@ jobs: strategy: matrix: os: [ubuntu-latest] - python-version: [3.6, 3.7, 3.8] + python-version: ['3.7', '3.8', '3.9', '3.10', '3.11'] + include: + - os: ubuntu-20.04 + python-version: '3.6' steps: - - uses: actions/checkout@v2 + - uses: actions/checkout@v3 - name: Set up Python ${{ matrix.python-version }} - uses: actions/setup-python@v2 + uses: actions/setup-python@v4 with: python-version: ${{ matrix.python-version }}