// 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 lambdaDaughter //daughter 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