<jupyter_start><jupyter_text>Some explanations What is the pair correlation function?In very short: the pair correlation function (aka radial distribution function) represents "the density of points found in average at a distance r of a given point in a sample".See the [Wikipedia page](https://en.wikipedia.org/wiki/Radial_distribution_function) for more information (but I assume that you already are familiar with this notion if you are looking for a script that computes it). What do you mean by "corrected to take account of boundary effects"?When the set of points is of finite size (which often occurs), some points near the boundaries will have less neighbours than what they would have in an infinite sample. This effect has to be corrected to properly compute the pair correlation function.This scripts has two methods to deal with boundaries:- the "normalization factors method": for the points that are too close to a boundary, the number of other points found at a given distance from them will be corrected to take account of the fact that in the bulk those particles should have more neighbours. This is the default behavior: all the points are considered (i.e. no data is lost), at the price of a time consuming computation.- the "exclusion method": all the points that are too close to a boundary are simply excluded from the computation. In this case, the computation is faster, at the price of dropping some points (i.e. losing data).  A simple example:<jupyter_code>import numpy as np
import matplotlib.pyplot as plt

from paircorrelation2d import pcf2d<jupyter_output><empty_output><jupyter_text>Create a hexagonal-like array of pointsThis will be the set of points for which we want to compute the pair correlation function in this example<jupyter_code>l_size=100 #the points will be placed in a square of size l_size*l_size
noise_amp=0.25 #we add some noise to mimic "real" data

col=np.arange(l_size)
points=np.zeros((l_size*l_size,2))

noise=np.random.rand(l_size*l_size,2)*noise_amp

for ii in range(l_size):
    points[ii*l_size:(ii+1)*l_size,0]=col+np.ones(l_size)*(1+(-1)**ii)/4+noise[ii*l_size:(ii+1)*l_size,0]
    points[ii*l_size:(ii+1)*l_size,1]=np.ones(l_size)*ii+noise[ii*l_size:(ii+1)*l_size,1]<jupyter_output><empty_output><jupyter_text>Let's look at the set of points:<jupyter_code>plt.scatter(points[:,0],points[:,1])<jupyter_output><empty_output><jupyter_text>Let's look at a subset of points (to see the general pattern):<jupyter_code>plt.scatter(points[:,0],points[:,1])
plt.axis([10,20,10,20])<jupyter_output><empty_output><jupyter_text>Compute the pair correlation function (pcf) taking account of all points:<jupyter_code>bins=np.linspace(0,5,100) #Since the distance between two particles is about 1 here, 
#we choose to compute the pcf only for distances up to 5.

[g_of_r_all,r] = pcf2d(points,bins,show_timing=True)
#the "show_timing" argument let you knows how long the script takes to run<jupyter_output>Creating boundary polygon and array of points inside took 0.063812 s
Creating all ring polygons took 0.015636 s
Computing normalization factors took 34.886039 s
Computing g(r) took 8.330064 s
Total time: 43.295551 s for 10000 points<jupyter_text>Here you can see that the total computation time is about 45 s (and that the computation of the normalization factor is the most time consuming operation).<jupyter_code>plt.plot(r,g_of_r_all)
plt.xlabel('r')
plt.ylabel('g(r)')
plt.title('Taking account of all points')<jupyter_output><empty_output><jupyter_text>Compute the pair correlation function (pcf) excluding points too close to the boundary:<jupyter_code>[g_of_r_exclude,r] = pcf2d(points,bins,fast_method=True,show_timing=True)
#the "show_timing" argument let you knows how long the script takes to run<jupyter_output>Creating boundary polygon and array of points inside took 0.063288 s
Creating all ring polygons took 0.021129 s
Computing normalization factors took 0.131086 s
Computing g(r) took 5.544017 s
Total time: 5.759521 s for 8090 points<jupyter_text>Here you can see that the total computation time is about 6 s (and that the computation of g(r) is faster than in the previous case because less points are considered).<jupyter_code>plt.plot(r,g_of_r_exclude)
plt.xlabel('r')
plt.ylabel('g(r)')
plt.title('Excluding points too close to the boundary')<jupyter_output><empty_output><jupyter_text>Wait, they look exactly the same, don't they?<jupyter_code>plt.plot(r,g_of_r_all,label='all points')
plt.plot(r,g_of_r_exclude,label='excluding points')
plt.xlabel('r')
plt.ylabel('g(r)')
plt.legend()<jupyter_output><empty_output><jupyter_text>This is because we have a set of points that is "bulky" (enough points are at a distance > 5 from the boundaries), but this is not always the case. Let's take another set of points:<jupyter_code>subpoints=points[np.where(points[:,1]<15)[0],:] #we exclude all the points with y>20

plt.scatter(subpoints[:,0],subpoints[:,1])
[g_of_r_all,r] = pcf2d(subpoints,bins)
[g_of_r_exclude,r] = pcf2d(subpoints,bins,fast_method=True)

plt.plot(r,g_of_r_all,label='all points')
plt.plot(r,g_of_r_exclude,label='excluding points')
plt.xlabel('r')
plt.ylabel('g(r)')
plt.legend()<jupyter_output><empty_output><jupyter_text>Here the difference is more noticeable, because more than half of the points are close enough to a boundary (i.e. their distance to the boundary is less than 5, which is the maximal distance we want for computing g(r)).  More information about boundaries:  You can define the boundary you want for your set of points:The scripts is based on the Polygon objects from [shapely](https://shapely.readthedocs.io/en/latest/manual.htmlpolygons) so any list of coordinates that creates a valid Polygon for shapely will work here. The scripts then automatically excludes the points that are not inside the area of interest you have defined.For example, if ones want to do a L-shape boundary:<jupyter_code>lshape_coord=np.array([[0,0],[100,0],[100,25],[50,25],[50,100],[0,100]])

#below is just for illustration purpose (it's not needed otherwise)
plt.plot(np.append(lshape_coord[:,0],0),np.append(lshape_coord[:,1],0))<jupyter_output><empty_output><jupyter_text>You can now compute the g(r) for the points of your set that are inside this L-shape polygon (and you can verify what points are kept by using the "plot=True" option):<jupyter_code>[g_of_r_lshape,r]=pcf2d(points,bins,coord_border=lshape_coord,plot=True)<jupyter_output><empty_output><jupyter_text>You can add holes in your area of interest:This might be useful for example if you are looking at a set of particles coordinates that are in a geometry with obstacles (the positions of the obstacles are exclusion zones where no particles can ever be found, so you have to remove them from the area of interest).<jupyter_code>square_coord=np.array([[0,0],[0,100],[100,100],[100,0]])
holes_coord=[np.array([[10,10],[10,30],[30,10]]),np.array([[60,60],[60,80],[80,80],[80,60]])] #the coordinates of holes polygons must be into a list (even when there's only one)

#below is just for illustration purpose (it's not needed otherwise)
plt.plot(np.append(square_coord[:,0],0),np.append(square_coord[:,1],0))
plt.plot(np.append(holes_coord[0][:,0],10),np.append(holes_coord[0][:,1],10),'r')
plt.plot(np.append(holes_coord[1][:,0],60),np.append(holes_coord[1][:,1],60),'r')<jupyter_output><empty_output><jupyter_text>You can now compute the g(r) for the points of your set that are inside this area of interest with holes (and you can verify what points are kept by using the "plot=True" option):<jupyter_code>[g_of_r_holes,r]=pcf2d(points,bins,coord_border=square_coord,coord_holes=holes_coord,plot=True)<jupyter_output><empty_output><jupyter_text>Two things to keep in mind about boundaries:- When no boundary is provided, the scripts computes the minimal convex polygon containing all the points in array_positions (the  convex hull). If the set of point has a non-convex boundary, the g(r) that will be computed will be wrong. For example, the convex hull of a L-shape set of points looks like this (polygon in blue, convex hull in red):<jupyter_code>plt.plot(np.append(lshape_coord[:,0],0),np.append(lshape_coord[:,1],0))
plt.plot([0,100,100,50,0,0],[0,0,25,100,100,0],'r')<jupyter_output><empty_output><jupyter_text>- The list of coordinates you provide for the boundary of the area of interest has to be "valid" in the sens used by the shapely library: linking all the points in order should result in a simple polygon with no line intersecting each other. For example:<jupyter_code>valid_square = np.array([[0,0],[0,100],[100,100],[100,0]])
plt.plot(np.append(valid_square[:,0],0),np.append(valid_square[:,1],0))<jupyter_output><empty_output><jupyter_text>This is a valid polygon.<jupyter_code>invalid_square = np.array([[0,0],[0,100],[100,0],[100,100]])
plt.plot(np.append(invalid_square[:,0],0),np.append(invalid_square[:,1],0),'r')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10839865621500559, 'fraction_numerical': 0.033370660694288914, 'mean_word_length': 6.889575971731449, 'pattern_counts': {'":': 2, '<': 68, '<?xml version=': 0, '>': 69, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '203853', 'license': 'MIT', 'n_tokens_mistral': 3114, 'n_tokens_neox': 2855, 'n_words': 982, 'path': 'example.ipynb', 'repo_name': 'aberut/paircorrelation2d', 'size': '8930'}