/*Information
This macro asks for a text file containing the data output by the psFRET_dual_view_image_analysis macro
It is currently expecting images collected using Micromanager simply to generate a BaseFileName hence "_MMStack" in one of the lines below.
Most of our data has three channels with the photoswitching data in Channel3.  This is why you'll see booleans checking for "Channel3"
*/

macro"PSFRET curve fit"{
cameraGain=2.17;
binFactor=4;
lambda=500;//wavelenth in nm
NA=1.4;
pixSize=39;//nm
pixCount=5000;
print("\\Clear");
run("Set Measurements...", "redirect=None decimal=9");
run("Clear Results");
OpenFilePath=File.openDialog("Pick the txt file with the data");
ImportResultsTable(OpenFilePath);
BaseFileName=substring(File.nameWithoutExtension,0,indexOf(File.nameWithoutExtension,"_MMStack"));
dir0=File.directory;
print(dir0);
print(BaseFileName);

headings=split(String.getResultsHeadings);
exptsToFit=newArray();
timesForFit=newArray();
for(i=2;i<headings.length;i++){
	if(findStringFragment(headings[i],"time")==false){
	exptsToFit=appendToArray(headings[i], exptsToFit);
	timesForFit=appendToArray(toString(headings[i]+"_time"), timesForFit);
	}
}


single_exp_fit="y = a*exp(-b*x)+c";

Dialog.create("Input cycle information");
Dialog.addMessage("*****************************************************************************");
Dialog.addNumber("Number_of_cycles",3);
Dialog.addNumber("Images_per_cycle",300);
Dialog.addMessage("");
Dialog.addCheckbox("Show plots and fits?", false);
Dialog.show();

numOfCycles=Dialog.getNumber();
numImagesPerCycle=Dialog.getNumber();
showPlots=Dialog.getCheckbox();

decayFitParametersArray=newArray("ExptCellNum"+","+"CycleNumber"+","+"A1"+","+"k"+","+"offset"+","+"t1/2"+","+"R2"+","+"Chi2");

for(i=0;i<exptsToFit.length;i++){
	expt=BaseFileName+"_"+exptsToFit[i];
	decayFitParametersArray=fitPSFRETDecayData(numOfCycles, numImagesPerCycle, expt, exptsToFit[i], timesForFit[i], decayFitParametersArray);
}

run("Clear Results");

exptDataArray=newArray("DDpre"+","+"DDpost"+","+"DAPre"+","+"DApost");

for(i=0;i<exptsToFit.length;i++){
	exptDataArray=getExptData(numOfCycles, numImagesPerCycle, dir0, exptDataArray);
}


single_exp_fit_data=File.open(dir0+"single_exp_fit_data.csv");
for(k=0;k<decayFitParametersArray.length;k++){
print(single_exp_fit_data,decayFitParametersArray[k]+",,"+exptDataArray[k]);
}
File.close(single_exp_fit_data);

if(File.exists(dir0+BaseFileName+"_single_exp_fit_data.csv")){
File.delete(dir0+BaseFileName+"_single_exp_fit_data.csv");
File.rename(dir0+"single_exp_fit_data.csv",dir0+BaseFileName+"_single_exp_fit_data.csv");
}else{
File.rename(dir0+"single_exp_fit_data.csv",dir0+BaseFileName+"_single_exp_fit_data.csv");
}






function getExptData(numCycles, imagesPerCycle, dir, arrayToReturn){
	dirList=getFileList(dir);
	for(i=0;i<dirList.length;i++){
		if(findStringFragment(dirList[i],"_green_raw_values")==true && findStringFragment(File.getParent(dir+File.separator+dirList[i]),"Channel3")==true && findStringFragment(dirList[i],"_single_exp_fit_")!=true)
		greenFilePath=dir+File.separator+dirList[i];
		if(findStringFragment(dirList[i],"_red_raw_values")==true && findStringFragment(File.getParent(dir+File.separator+dirList[i]),"Channel3")==true && findStringFragment(dirList[i],"_single_exp_fit_")!=true)
		redFilePath=dir+File.separator+dirList[i];
	}
	ImportResultsTable(greenFilePath);
	greenExptData=parseResultsTable(numCycles, imagesPerCycle);
	run("Clear Results");

	ImportResultsTable(redFilePath);
	redExptData=parseResultsTable(numCycles, imagesPerCycle);

	selectWindow("Results");
	run("Close");
	for(i=0;i<greenExptData.length;i++){
		arrayToReturn=appendToArray(toString(greenExptData[i])+","+toString(redExptData[i]),arrayToReturn);
	}

return arrayToReturn;
}

function parseResultsTable(numCycles, imagesPerCycle){
	headings=split(String.getResultsHeadings);
	exptDataToFind=newArray();
	for(i=2;i<headings.length;i++){
	/*for(i=0;i<headings.length;i++){*/ //use this when no background image
		if(findStringFragment(headings[i],"time")==false){
		exptDataToFind=appendToArray(headings[i], exptDataToFind);
		}
	}
	exptDataToReturn=newArray();
	
	for(i=0;i<exptDataToFind.length;i++){
		for(j=0;j<numCycles;j++){		
			PrePS=getResult(exptDataToFind[i],(j*imagesPerCycle));
			//print(PrePS);
			PostPS=getResult(exptDataToFind[i],(j*imagesPerCycle)+(imagesPerCycle-1));
			//print(PostPS);
			exptDataToReturn=appendToArray(toString(PrePS+","+PostPS),exptDataToReturn);			
		}
	}
return exptDataToReturn;
print(exptDataToReturn);
}

function fitPSFRETDecayData(numCycles, imagesPerCycle, exptN, exptArray, timeArray, arrayToReturn){

	for(i=0;i<numCycles;i++){
		x=newArray(imagesPerCycle);
		y=newArray(imagesPerCycle);		
		for(j=0;j<x.length;j++){
			x[j]=getResult(timeArray,(i*imagesPerCycle)+j);
			y[j]=getResult(exptArray,(i*imagesPerCycle)+j);
		}
		Array.getStatistics(x,xminVal,xmaxVal,xmeanVal,xsdVal);
		x=arraySubtractValue(x,xminVal);
		single_exp_fit="y = a*exp(-b*x)+c";
		t12=findT1half(x,y);
		if(t12==isNaN(parseFloat(t12)))
			t12=2;
		initialGuesses=newArray(y[0],1/t12,y[y.length-1]);
		Fit.doFit(single_exp_fit,x,y,initialGuesses);
		a=Fit.p(0);
		b=Fit.p(1);
		c=Fit.p(2);
		R2=Fit.rSquared;
		theFit=getTheFit(x,a,b,c);
		Chi2=getChi2(y,theFit);
		if(showPlots){
			Plot.create("Data and fit","Time(sec)","Fluorescence");
			Plot.add("circles",x,y);
			Plot.setColor("red");
			Plot.add("line",x,theFit);
			Plot.addText("a="+a,0.1,0.1);
			Plot.addText("k="+b,0.1,0.15);
			Plot.addText("offset="+c,0.1,0.2);
			Plot.addText("Chi2="+Chi2,0.1,0.25);
			Plot.show();
		}	
		fitParameters=toString(exptN+","+i+","+a+","+b+","+c+","+t12+","+R2+","+Chi2);
		arrayToReturn=appendToArray(fitParameters,arrayToReturn);
	}
return arrayToReturn;
}//end function fitPSFRETDecayData

function getTheFit(time,A,k,C){
	arrayToReturn=newArray(time.length);
	for(i=0;i<time.length;i++){
		arrayToReturn[i]=(A*exp(-k*time[i]))+C;
	}
	return arrayToReturn;
}

function getChi2(data, fit){
	AiryRadius = 1.22*lambda/NA;         
    PSFArea = PI*AiryRadius*AiryRadius;        
    numPixelsInPSF = PSFArea/(pixSize*pixSize);
	s2=0;     
	sumWeightedVariance=0;
	divideArrayByValue(fit,cameraGain);
	divideArrayByValue(data,cameraGain);
	data1=newArray(data.length);
	for(i=0;i<data.length;i++){
		data1[i] = data[i];
	}

	multiplyArrayByValue(fit,pixCount);
	multiplyArrayByValue(data1,pixCount);
	multiplyArrayByValue(data,pixCount*binFactor*binFactor);
	coVarArray=newArray(data.length);
	for(i=0;i<data.length;i++){
		coVarArray[i] = data[i]*(pixCount*binFactor*binFactor/numPixelsInPSF);
	}
	for(i=0;i<data.length;i++){
		data[i] = data[i]+coVarArray[i];
	}
	            
	for(i=0;i<data.length;i++){
		sumWeightedVariance += (1/data[i]);
	}
	averageWeightedVariance=sumWeightedVariance/data.length;
	weightArray = newArray(data.length);
	
	for(i=0;i<data.length;i++){
		weightArray[i] += (1/data[i])/averageWeightedVariance;
	}
	for(ii=0;ii<data.length;ii++){
		s2 += ((pow((data1[ii]-fit[ii]),2))*weightArray[ii]);
	}
	s2=s2/(data.length-4);
	Chi2ToReturn=s2*averageWeightedVariance;
	return Chi2ToReturn;
}

 function ImportResultsTable(filePath) {
     requires("1.35r");
     lineseparator = "\n";
     cellseparator = ",\t";

     // copies the whole RT to an array of lines
     lines=split(File.openAsString(filePath), lineseparator);

     // recreates the columns headers
     labels=split(lines[0], cellseparator);
     if (labels[0]==" ")
        k=1; // it is an ImageJ Results table, skip first column
     else
        k=0; // it is not a Results table, load all columns
     for (j=k; j<labels.length; j++)
        setResult(labels[j],0,0);

     // dispatches the data into the new RT
     run("Clear Results");
     for (i=1; i<lines.length; i++) {
        items=split(lines[i], cellseparator);
        for (j=k; j<items.length; j++)
           setResult(labels[j],i-1,items[j]);
     }
     updateResults();
 }

function findT1half(xArrayForT1half,yArrayForT1half){
Array.getStatistics(xArrayForT1half,xmin,xmax,xmean,xsd);
Array.getStatistics(yArrayForT1half,ymin,ymax,ymean,ysd);
t1halfToReturn=0;
i=0;
while(yArrayForT1half[i]>ymax/2){
	i++;
	if(i==yArrayForT1half.length)
	return "noTOneHalf";
}
t1halfToReturn=xArrayForT1half[i]-xmin;
return t1halfToReturn;
}

function arraySubtractValue(theArray,value){
for(i=0;i<theArray.length;i++){
theArray[i]=theArray[i]-value;
}
return theArray;
}

function divideArrayByValue(theArray,value){
for(i=0;i<theArray.length;i++){
theArray[i]=theArray[i]/value;
}
return theArray;
}

function multiplyArrayByValue(theArray,value){
for(i=0;i<theArray.length;i++){
theArray[i]=theArray[i]*value;
}
return theArray;
}

function findStringFragment(theString,theFragment){
	stringFound=false;
for(i=0;i<=lengthOf(theString)-lengthOf(theFragment);i++){
	check4Frag=substring(theString,i,i+lengthOf(theFragment));
	if(check4Frag==theFragment)
	stringFound=true;
}
	return stringFound;
}

function appendToArray(value, array) {
temparray=newArray(lengthOf(array)+1);
for(i=0;i<lengthOf(array);i++){
temparray[i]=array[i];
}
temparray[lengthOf(temparray)-1]=value;
array=temparray;
return array;
}

}//end of macro