##### https://github.com/GrandlLab/PiezoClusteringModel
Tip revision: c446386
Piezo Clustering Model
``````// this code was inspired by “Simulating a Thomas cluster point process”  in R from Paul Keeler;
// https://hpaulkeeler.com/simulating-a-thomas-cluster-point-process/
// https://github.com/hpaulkeeler/posts/tree/master/ThomasClusterRectangle

This was written and tested in Igor Pro 8.0 (Wavemetrics)

#pragma TextEncoding = "UTF-8"
#pragma rtGlobals=3		// Use modern global access method and strict wave access.

function optimizeclustering(lambdaparent,lambdadaughter,sigma,numb_samplepoints,datainput) //this function creates simulated points via a Thomas point process, calculates histogramss for simulated and real data, and calculates the sum of squares for the difference between the two.

variable lambdaParent //parent points
variable sigma  //sigma
variable numb_samplepoints/ /number of sample centers- each is the center of a "patch" with a given radius.
string datainput  //this is the raw data wave. Current amplitudes in A, one value per cell

variable xmin=0,xmax=50,ymin=0,ymax=50 //setting up the boundaries of the big box in microns. The full box will be twice the width (here, -50 to +50).
variable xdelta=xmax-xmin, ydelta=ymax-ymin //calculating the length of the box sides
variable areatotal = xdelta*ydelta //calculatingthe area of the box, for setting up density of points

variable rExt=7*sigma //This accounts for edge effects by extending the original box slightly, based on value of sigma, so all points end up in box.
variable xMinExt=xMin-rExt, xMaxExt=xMax+rExt, yMinExt=yMin-rExt, yMaxExt=yMax+rExt //acutal values for extended box

variable xDeltaExt=xMaxExt-xMinExt, yDeltaExt=yMaxExt-yMinExt //extended box sides
variable areaTotalExt=xDeltaExt*yDeltaExt //extended box area

variable numbPointsParent = (areaTotalExt*lambdaparent) //calculates number of parent points- based on density and area.

Make/O/N=(numbPointsParent) xxparent=enoise(xdeltaext) //assigns x coordinates to parent points from a uniform distribution
Make/O/N=(numbPointsParent) yyparent=enoise(ydeltaext) //assigns y coordinates to parent points from a uniform distribution
Make/O/N=(numbPointsParent)  numbPointsDaughter = poissonnoise(lambdaDaughter) //assigns a # of daughter points to each parent point, poisson distributed. Note that to get a truly random distribution, simply set numbpointsdaughter=1
variable  numbPoints=sum(numbPointsDaughter) //calculates total number of daughter points (Sum of each assigned to a given parent)

Make/O/N=(numbPoints) xx0=gnoise(sigma) //essentially a "noise" wave, the length of the total number of daughter points, with stdev of sigma

Make/O/N=(numbPoints) yy0=gnoise(sigma) //same as above

//this whole next section makes the final daughter coordinates by first replicating the xx coordinates of each parent point n times, where n is the number of daughter points assigned to that parent.
Make/O/N=(numbPoints) xx=0, yy=0, xxnoiseless=0, yynoiseless=0
variable xxwaveindex = 0
xxnoiseless[0,(sum(numbPointsDaughter,0,0))]=xxparent[0]
yynoiseless[0,(sum(numbPointsDaughter,0,0))]=yyparent[0]
xxwaveindex+=1
variable startpoint
variable endpoint

do
startpoint=sum(numbpointsdaughter,0,(xxwaveindex-1))
endpoint = sum(numbpointsdaughter,0,xxwaveindex)
xxnoiseless[startpoint,endpoint] = xxparent[xxwaveindex]
yynoiseless[startpoint,endpoint] = yyparent[xxwaveindex]
xxwaveindex += 1
while (xxwaveindex <= (numbPointsparent))

xx=xxnoiseless+xx0 //adding in noise in xx coordinates
yy=yynoiseless+yy0 //adding in noise in yy coordinates

concatenate/O {xx,yy}, simchannelmatrix //putting the xx and yy coordinates together in a matrix.

variable x2min = (-xmax) //setting up bounds for sample point box.
variable x2max = (xmax)
variable y2min = (-ymax)
variable y2max = (ymax)
variable x2delta=x2max-x2min
variable y2delta=y2max-y2min
variable area2total = x2Delta*y2Delta

Make/O/N=(numb_samplepoints) xx_sample=enoise(.5*x2delta) //xx locations of pipette centers from a uniform distribution
Make/O/N=(numb_samplepoints) yy_sample=enoise(.5*y2delta) //yy locations of pipette centers from a uniform distribution

concatenate/O {xx_sample,yy_sample}, samplematrix //now the sample points are in a matrix.

Matrixop/o/free matAx=colRepeat(xx,numb_samplepoints)//matrix of channel x coordinates
Matrixop/o/free matAy=colRepeat(yy,numb_samplepoints)//matrix of channel y coordinates
Matrixop/o/free matBx=rowRepeat(xx_sample,numbpoints) //matrix of pipette center x coordinates
Matrixop/o/free matBy=rowRepeat(yy_sample,numbpoints)//marix of pipette center y coordinates

MatrixOP/o distances=sqrt((matAx-matBx)*(matAx-matBx)+(matAy-matBy)*(matAy-matBy)) //distance matrix between all channel points and all pipette points.

//Now making a matrix of normally distributed pipette sizes
Make/O/N=(numb_samplepoints) pipettesize=((0.8+gnoise(.147))) //average pipette radius is 0.8 um, with a standard deviation of .147 um
Matrixop/o pipettematrix=rowrepeat(pipettesize,numbpoints)

//now counting channels that fall within the pipette radii.
MatrixOP/O clusterdistances=greater(pipettematrix,distances) //compares distance matrix to pipette matrix and assigns a binary value (0 for distance > pipette radius, 1 for distance < pipette radius)
MatrixOP/O/NPRM/S channelcount= sumcols(clusterdistances) //sums the columns of above matrix- i.e., # of distances that are within the cutoff of pipette radius
MatrixOP/O channelcount = redimension(channelcount,numb_samplepoints,1) //redimensions above matrix into a column for further analysis

//converting channel counts to currents
Make/O/R/N=(numb_samplepoints) simchannelampnoise= gnoise(2.22e-13)+9.76e-13 //this makes a normal distribution of single channel conductances based on N2A values- stdev of 2.45e-13 and mean 9.38e-13
MAke/O/R/N=(numb_samplepoints) simchannelamp = channelcount*simchannelampnoise //converts simulated channel counts to single channel conductances

duplicate/O \$(datainput) realdata //placeholder wave for raw data.
Make/O/N=100/O simchannelamp_Hist,realN2Acount //making two empty histograms for data comparison
Histogram/B={-.49e-12,.976e-12,100} simchannelamp,simchannelamp_Hist //binning simulated channel amplitudes- bin width is mean single channel conductance and offset by half that value, so should be centered around 1 channel
Histogram/B={-.49e-12,.976e-12,100} realdata, realN2Acount //binning real chanel amplitudes
duplicate/O simchannelamp_hist simchannelnorm //creating normalized wave for simulated channel counts
duplicate/O realN2Acount realdatanorm //creating normalized wave for real data
simchannelnorm=simchannelnorm/sum(simchannelamp_hist)  //calculating normalized wave
realdatanorm = realn2acount/sum(realn2acount) //calculating normalized wave
Make/O/N=100 histogramdiff //making error wave
histogramdiff = (simchannelnorm-realdatanorm)^2 //calculating square difference for each bin
variable sqdiff = sum(histogramdiff) //sum of squared differences

return sqdiff // IMPORTANT: this "returns" this value from this function- so if you call this function from another function for fitting purposes, this is the output.

end

``````