(* This script scans through the actin distributions for a set of experiments and for each distribution (1) locates the peak intensity in the raw intensity, (2) locates the peak moment in a smoothed intensity, (3) calculates the second moments about the peak in the smoothed intensity, and (4) and calculates the Fisher information (in both the raw distribution) *)

(* Set up string variables needed to read the distributions. 
NOTE: this may need to be adjusted by the user to account for different file naming conventions *)

neuronExperimentTag = "ABL KD"  ; timeTag ="3min" ; numberOfNeurons = 14 ;

(* neuronExperimentTag = "ABL KD"  ; timeTag ="1min" ; numberOfNeurons = 13 ; *)

(* neuronExperimentTag = "ABL GOF" ;  timeTag ="3min" ; numberOfNeurons = 14 ; *)

(* neuronExperimentTag = "ABL GOF" ;  timeTag ="1min" ; numberOfNeurons = 13 ; *)

(* neuronExperimentTag = "WT" ; timeTag ="3min" ; numberOfNeurons = 14 ; *)

(* neuronExperimentTag = "WT" ; timeTag ="1min" ; numberOfNeurons = 14 ; *)

(* neuronExperimentTag = "Dicer" ; timeTag ="3min" ; numberOfNeurons = 14 ; *)

neuronNameList =  Table[   neuronExperimentTag  <> " " <> ToString[ nx ], { nx, 1, numberOfNeurons } ]

(* Set up some paths *)

pathToTraceDirectory =  "/Users/mcqueenp/ActiveProjects/AxonGrowthInDrosophilaWing/TraceData/" <> If[ neuronExperimentTag == "WT", "Wild Type",  neuronExperimentTag  ] <> " " <> timeTag <> "/"

pathToOutputDirectory =  "/Users/mcqueenp/Downloads/"

timeInfoFile = "IntervalBetweenFrames.csv"  ;

beginningPointOnAxonFile = "BeginningMicrometer.csv"   ;

actinFileSuffix= "_actin.csv"   ;

verySmallNumber = 10^-10

binWidth = 0.1

windowSize = 5.0

binsInHalfWindow = Ceiling[ windowSize/( 2*binWidth ) ]

timeInfoFile = "IntervalBetweenFrames.csv"

beginningPointOnAxonFile = "BeginningMicrometer.csv"

actinFileName = "actin.csv"

(* Set path to directory with individual neuron data *)

neuronNameAll = {} ; timeIndexAll = {} ; timeSinceBeginningAll = {} ; positionMaxRawAll = {} ;   fisherInformationAll = {} ;  positionMaxSmAll = {} ; secondMomentSmAll = {} ;secondMomentBodySmAll = {} ; secondMomentEndSmAll = {} ; filesWithFailedImports = {} ;

For[ nx = 1, nx ≤ Length[ neuronNameList ], nx++,
pathToNeuronDirectory = pathToTraceDirectory <> neuronNameList[[ nx ]] <> "/" ;
(* Import time information data *)
fileName = pathToNeuronDirectory <> timeInfoFile ;
fileThere = FileExistsQ[ fileName ] ; 
If[ !fileThere, AppendTo[ filesWithFailedImports, fileName ]  ] ;
If[ !fileThere, Continue[ ] ] ;
dumArray = Import[ fileName ] ;
maximumTimeIndex = dumArray[[ 2 ]][[ 2 ]] ;
minimumTimeIndex = dumArray[[ 2 ]][[ 1 ]] ; deltaTime = ( 60.0*dumArray[[ 2 ]][[ 3 ]] + dumArray[[ 2 ]][[ 1 ]])/60.0 ;
dumArray =. ;


fileName = pathToNeuronDirectory <>  beginningPointOnAxonFile ;
fileThere = FileExistsQ[ fileName ] ;
If[ fileThere, dumArray = Import[ fileName ], AppendTo[ filesWithFailedImports, fileName ]   ] ;
If[ fileThere,beginningPointForAnalysis = dumArray[[ 2 ]][[ 1 ]], beginningPointForAnalysis = 0.0 ] ;
If[ fileThere, dumArray =. ] ;




(* Import the actin data for each time point *)
For[ tx = minimumTimeIndex , tx ≤ maximumTimeIndex,  tx++,


timeString ="t=" <> ToString[ tx ] <> "_"  ;

fileName = pathToNeuronDirectory <> timeString <> actinFileName  ;
fileThere = FileExistsQ[ fileName ] ;
If[ !fileThere, AppendTo[ filesWithFailedImports, fileName ]  ] ;
If[ !fileThere, Continue[ ] ] ;


dumArray= Import[ fileName ] ;
If[ dumArray == $Failed, AppendTo[ filesWithFailedImports, fileName ]  ] ;
If[ dumArray == $Failed, Continue[ ] ] ;
(* Skip this time point if the data for it cannot be found *)

actinData= SortBy[ Table[ { dumArray[[ n ]][[ 4 ]], dumArray[[ n ]][[ 5 ]] }, { n, 5, Length[ dumArray ] } ], First ] ;

AppendTo[ neuronNameAll, neuronNameList[[ nx ]]  ] ;
AppendTo[ timeIndexAll, tx ] ;

AppendTo[ timeSinceBeginningAll,
If[ minimumTimeIndex == 0, deltaTime* tx, deltaTime*( tx - 1 ) ] ] ;

dumArray =. ;

edgeLength= Max[ Table[  actinData[[ ix  ]][[ 1 ]], { ix, 1,  Length[ actinData ] }  ]  ] ;

(* Find the location of the maximum actin (raw) intensity for this time point *)
maxactRaw = -1000000.0 ; posMaxAct = -1000000.0 ;
For[ ix = 1, ix ≤ Length[ actinData  ], ix++,

If[ beginningPointForAnalysis  > actinData[[ ix ]][[ 1 ]], Continue[ ] ] ;

If[ maxactRaw > actinData[[ ix ]][[ 2 ]] , Continue[ ] ] ;
  maxactRaw = actinData[[ ix ]][[ 2 ]] ;
 posMaxAct = actinData[[ ix ]][[ 1 ]] ;

] ;

AppendTo[ positionMaxRawAll, posMaxAct ] ;


(* Find the trapezoid approximation to the actin integral. This is needed for determining the Fisher Information *)
actinIntegralTrapApprx =  0.0 ; 
For[ ix = 1, ix < Length[ actinData ], ix++, 
If[ beginningPointForAnalysis  > actinData[[ ix ]][[ 1 ]], Continue[ ] ] ;
deltaX = actinData[[ ix + 1 ]][[ 1 ]] - actinData[[ ix ]][[ 1 ]] ; 
actLocal =  actinData[[ ix + 1 ]][[ 2 ]] + actinData[[ ix ]][[ 2 ]] ; 
actinIntegralTrapApprx  += 0.5*deltaX*actLocal ;
] ;



(* Find the Fisher information in the raw distribution *)
uuu = {} ; 
For[ ix = 1, ix ≤ Length[ actinData ], ix++, 
AppendTo[ uuu, If[ actinData[[ ix ]][[ 2 ]] > 0.0, 
( actinData[[ ix ]][[ 2 ]]/actinIntegralTrapApprx )^0.5, 
0.0 ] 
] ;
] ; 

 xxx = {}  ;
derivUUU = {}  ;
For[ ix =  2, ix ≤ Length[ uuu ], ix++, 

If[ beginningPointForAnalysis  > actinData[[ ix - 1 ]][[ 1 ]], Continue[ ] ] ;

derv = If[ actinData[[ ix ]][[ 1 ]] == actinData[[ ix - 1 ]][[ 1 ]], 0, ( uuu[[ ix ]] - uuu[[ ix - 1 ]] )/( actinData[[ ix ]][[ 1 ]] - actinData[[ ix - 1 ]][[ 1 ]] ) ] ;
AppendTo[ derivUUU, derv ] ;
AppendTo[ xxx, 0.5*(  actinData[[ ix ]][[ 1 ]] + actinData[[ ix - 1 ]][[ 1 ]] ) ] ;
 ] ;

fInf = 0.0 ;
For[ ix = 2, ix ≤ Length[ xxx ], ix++, 

deltaX  =  xxx[[ ix ]] - xxx[[ ix - 1 ]]  ;

qqqA =derivUUU[[ ix - 1 ]]^2  ;
qqqB = derivUUU[[ ix ]]^2 ; 
fInf += 0.5*deltaX*( qqqA + qqqB ) ;

] ; (* End loop over ix *)

fInf *= 4.0 ;
AppendTo[ fisherInformationAll, fInf ] ;

xxx=.;
uuu =. ;
derivUUU =. ;




(* Find the maximums in actin distributions smoothed over 5mum windows. Also find the second moments about that maximum *)

(* Create a binned distribution for the actin. First, make binned version of raw distribution *)
numberOfBins = 1 + Floor[ ( edgeLength- beginningPointForAnalysis )/binWidth ] ;

binAux = Table[ 0, { bx, 1, numberOfBins } ] ;

For[ inx = 1, inx ≤ Length[ actinData  ] - 1, inx++, 

If[  beginningPointForAnalysis > actinData[[ inx ]][[ 1 ]], Continue[ ] ]   ; 

targetA = actinData[[ inx ]][[ 1 ]] - beginningPointForAnalysis ;
binA = 1 + Floor[ targetA/binWidth ] ;
If[ binA < 1, binA = 1 ] ;
If[ binA > numberOfBins, binA = numberOfBins ] ;

binAux[[ binA ]] += actinData[[ inx ]][[ 2 ]] ;

] ;

(* Now smooth the data *)
actinSmoothed = Table[ 0,{ bx, 1, numberOfBins } ] ;

For[ bx = 1, bx ≤  numberOfBins, bx++, upperL = bx + binsInHalfWindow ;
If[ upperL > numberOfBins, upperL = numberOfBins ] ;
lowerL = bx - binsInHalfWindow ;
If[ lowerL < 1, lowerL = 1 ] ;
sumSum = Sum[ binAux[[ ix ]], { ix, lowerL, upperL } ] ;
actinSmoothed[[ bx ]] =1.0*sumSum/( upperL - lowerL ) ;
] ;


maxInSm = -1000000.0 ; indexOfMax = -1000000 ;
For[ inx = 1, inx ≤ Length[ actinSmoothed ], inx++, 
If[ maxInSm > actinSmoothed[[ inx ]], Continue[ ] ] ;
 
maxInSm = actinSmoothed[[ inx ]] ;
indexOfMax = inx ; 

] ;

AppendTo[ positionMaxSmAll, binWidth*indexOfMax + beginningPointForAnalysis ] ;


(* Find second moments *)
massSum = 0.0 ;  bodyMomentSquared = 0.0 ; endMomentSquared = 0.0 ;
For[ bx = 1, bx ≤ numberOfBins, bx++, 

massSum += actinSmoothed[[ bx ]] ;
If[ bx < indexOfMax, bodyMomentSquared += actinSmoothed[[ bx ]]*1.0*( indexOfMax - bx )^2 ] ;
If[ bx > indexOfMax, endMomentSquared += actinSmoothed[[ bx ]]*1.0*( bx - indexOfMax  )^2 ] ; 

] ;


bodyMomentSquared  = If[ massSum  > 0.0, bodyMomentSquared/massSum, 0.0 ] ;
endMomentSquared = If[ massSum > 0.0, endMomentSquared/massSum, 0.0 ] ;
secondMomentSquared = bodyMomentSquared +endMomentSquared ;

AppendTo[ secondMomentBodySmAll, binWidth*( bodyMomentSquared )^0.5 ] ;
AppendTo[ secondMomentEndSmAll, binWidth*( endMomentSquared )^0.5 ] ;
AppendTo[ secondMomentSmAll, binWidth*( secondMomentSquared )^0.5 ] ;

binAux =. ;
actinSmoothed =.;

actinData =.;

] ; 

]   (* End of loop on nx *)

(* Check on size of arrays. They all should be equal. *)

{ Length[ neuronNameAll ],   Length[timeIndexAll ],   Length[ timeSinceBeginningAll ],   Length[positionMaxRawAll ],     Length[ fisherInformationAll ],   Length[  positionMaxSmAll ],   Length[ secondMomentSmAll ],  Length[ secondMomentBodySmAll ],   Length[ secondMomentEndSmAll ] }

(* Which files did not get open, if any? *)

filesWithFailedImports

(* Export the computed results into a CSV file *)

exportArray = Table[ { neuronNameAll[[ elx ]], timeIndexAll[[ elx ]] ,timeSinceBeginningAll[[ elx ]],  positionMaxRawAll[[ elx ]], fisherInformationAll[[ elx ]], positionMaxSmAll[[ elx ]], secondMomentSmAll[[ elx ]], secondMomentBodySmAll[[ elx ]], secondMomentEndSmAll[[ elx ]]}, { elx, 1,Length[ neuronNameAll ] } ] ;

Export[  pathToOutputDirectory <>  neuronExperimentTag <> "_" <> timeTag <> "_Actin_Distribution_Properties.csv" ,  
Join[ { {  "Experiment Label", "Time Index", "Time since Beginning of experiment (sec)", "Position of maximum in raw actin intensity (mum)", "Fisher Information of raw intensity (1/mum)", "Position of maximum in smoothed intensity (mum)", "2nd Moment about maximum (smoothed intensity) (mum)", "Partial 2nd Moment towards body (mum)", "Partial 2nd Moment towards end (mum)"  } }, 
exportArray ], 
 "CSV" ]   ;

(* End of file *)