(*
	Example notebook for ImagingAnalysis.m (Rel. 1.0)

	2010/07/29
	Jihwan Myung
*)

(* Load ImagingAnalysis.m package *)

<<ImagingAnalysis`

(* Put ImagingAnalysis.m in (User Directory)/LIbrary/Mathematica/Applications *)

(* File import: Specify directories of data files. In the current example, preprocessed image shots *)

$ID="SampleDataFolder";

$HomeDir="SubFolder/"<>$ID<>"/";
ImageDirectory="HomeFolder/"<>$HomeDir<>"BmalLD-processed/";
RawDirectory="HomeFolder/"<>$HomeDir<>"BmalLD/";

(* CreateFullID makes database for each fileset  *)

id[1]=CreateFullID[$ID,ImageDirectory,RawDirectory,700,"LD"];

id[1][[1]]

(* Preprocessing *)

(* Onigiri sectioning defines basic tissue geometry. It works when there is a single slice with the dorsal side up in the viewfield *)

osec=OnigiriSection[id[1]];

(* sall is a detrended trace *)

sall=GetAllS[id[1]];

phaseds=PhaseDS[sall,id[1]];

rlall=GetAllRL[sall,id[1]];

rfall=GetAllRF[sall,id[1]];

FS=PhaseS[sall,id[1]];

abrac=Unwrap[FS];

abraclen=Table[Length[abrac[[n]]],{n,1,id[1][[7]]}];
exclude=Flatten[Position[abraclen,
	Complement[abraclen,{Median[abraclen]}][[1]]]];

mabrac=Mean[Table[abrac[[n]],{n,Complement[Range[1,id[1][[7]]],exclude]}]];
mphaseds=Mean[Table[phaseds[[n]],{n,Complement[Range[1,id[1][[7]]],exclude]}]];

Table[abrac[[exclude[[n]]]]=mabrac,{n,1,Length[exclude]}];
Table[phaseds[[exclude[[n]]]]=mphaseds,{n,1,Length[exclude]}];

PlotAllRawTimeSeries[id[1],SeriesTimeScale->Automatic,TimeLabel->"ExT",Backdrop->True]

Show[PlotAllTimeSeries[sall,id[1],SeriesTimeScale->Automatic,TimeLabel->"ExT",Backdrop->True],PlotRange->{{0,24*(60/id[1][[1]])+id[1][[2]]},{-Max[sall[[1]]],Max[sall[[1]]]}},AspectRatio->0.5]

96*4

(* sorting by mean over time *)

PlotSortedLuminescence[id[1],CalibrationBarLabel->"Raw luminescence",TimeLabel->"ExT"]

rawall=RawTimeSeries[id[1]];

fall=id[1][[9]];

sortindex=id[1][[8]];

SBright=Table[Total[Abs[sall[[n]]]],{n,1,id[1][[7]]}];

RBright=Table[Total[Abs[rawall[[n]]]],{n,1,id[1][[7]]}];

sbcl=FindClusters[RBright->Range[Length[RBright]]];

sbcl2=FindClusters[SBright->Range[Length[SBright]]];

sbcllabel=ClusterLabel[sbcl][[All,2]];

sbcllabel2=ClusterLabel[sbcl2][[All,2]];

CalibClusterPlot[sbcl,id[1],"Raw lum. cluster"]

CalibClusterPlot[sbcl2,id[1],"Luminescence cluster"]

fftper=Table[FFTPeriod[sall[[n]],id[1]],{n,1,id[1][[7]]}];

percl=FindClusters[fftper->Range[Length[fftper]]];

CalibClusterPlot[percl,id[1],"Period cluster"]

phadscl=FindClusters[phaseds->Range[Length[phaseds]]];

CalibClusterPlot[phadscl,id[1],"PhaseDS cluster"]

kin=KuramotoIndex[abrac,phadscl];

PlotPhaseCoherence[kin,id[1],phadscl]

kindex=KuramotoIndex[abrac];

PlotPhaseCoherence[kindex,id[1]]

Show[Table[ListPlot[Mean[phaseds[[phadscl[[n]]]]],PlotStyle->Hue[((n-1)/Length[phadscl]),1,0.8]],{n,1,Length[phadscl]}]]

pha=Table[Mean[phaseds[[n]]],{n,1,id[1][[7]]}];

phacl=FindClusters[pha->Range[Length[pha]]];

CalibClusterPlot[phacl,id[1],"Phase cluster"]

ratio=Table[Sum[If[rfall[[n]][[k]]≥0,rfall[[n]][[k]],0],{k,1,Length[rfall[[n]]]}]/
Sum[Abs[If[rfall[[n]][[k]]<0,rfall[[n]][[k]],0]],{k,1,Length[rfall[[n]]]}],{n,1,id[1][[7]]}];

ratiocl=FindClusters[ratio->Range[Length[ratio]]];

CalibClusterPlot[ratiocl,id[1],"Ratio cluster"]

ListPlot[rfall[[1]]]

corrmatrix=CorrMatrix[rfall,id[1]];

Zp=Z[FS];

pcorrmatrix=PartialCorrelation[rfall,Zp,corrmatrix];

nclust=NCluster[pcorrmatrix,id[1]]

ptbl=CorrCluster[pcorrmatrix,phaseds];

CalibClusterPlot[ptbl,id[1],"Corr. cluster"]

ntbl=CorrCluster[pcorrmatrix,phaseds,nclust];

CalibClusterPlot[ntbl,id[1],"Corr. cluster"]

Show[Table[ListPlot[Mean[phaseds[[ptbl[[n]]]]],PlotStyle->Hue[((n-1)/Length[ptbl]),1,0.8]],{n,1,Length[ptbl]}]]

kinp=KuramotoIndex[abrac,ptbl];

Show[PlotPhaseCoherence[kinp,id[1],ptbl],AspectRatio->0.5,
PlotRange->{{0,24*(60/id[1][[1]])+id[1][[2]]},{0,1.05}}]

ppos=PolarPosition[ptbl,osec,id[1]]

PlotPolarPosition[ppos,ptbl]

thrcl=FindClusters[Thread[{pha,fftper,ratio}]->Range[Length[ratio]]];

CalibClusterPlot[thrcl,id[1],"Complex cluster"]

ListPointPlot3D[Thread[{pha,fftper,ratio}]]

CLsall=Flatten[Table[sall[[ptbl[[n]]]],{n,1,Length[ptbl]}],1];

RasterPlot[CLsall,id[1],ptbl,CalibrationBarLabel->"Detrended luminescence",TimeLabel->"ExT",
MainLabel->"Ranked by cluster", Clustered->True]

CLrlall=Flatten[Table[rlall[[ptbl[[n]]]],{n,1,Length[ptbl]}],1];

RasterPlot[CLrlall,id[1],ptbl,CalibrationBarLabel->"Normalized luminescence",TimeLabel->"ExT",
MainLabel->"Ranked by cluster", Clustered->True]

CLrfall=Flatten[Table[rfall[[ptbl[[n]]]],{n,1,Length[ptbl]}],1];

RasterPlot[CLrfall,id[1],ptbl,CalibrationBarLabel->"Normalized luminescence",TimeLabel->"ExT",
MainLabel->"Ranked by cluster",Clustered->True]

pall=Table[Peaks[sall[[n]]],{n,1,id[1][[7]]}];

RasterPlotCounts[pall,id[1],Backdrop->True,TimeLabel->"ExT", Height->300,Clustered->False]

CLpall=Flatten[Table[pall[[ptbl[[n]]]],{n,1,Length[ptbl]}],1];

RasterPlotCounts[CLpall,id[1],ptbl,Backdrop->True,TimeLabel->"ExT", Height->300,Clustered->True]

sorr=ClusteredCorrMatrix[pcorrmatrix,ptbl];

ClusteredMatrixPlot[sorr,ptbl,"Partial correlation coefficient"]

ol=OnigiriLabel[osec];

CalibClusterPlot[osec[[Range[1,3]]],id[1],"Onigiri section"]

CalibRImagePlot[psl=Table[Mean[(12/Pi)phaseds[[n]]],{n,1,id[1][[7]]}],id[1],"Mean phasse dev. (hr)"]

CalibRImagePlot[fftper,id[1],"Period (hr)"]

CalibRImagePlot[psl=Table[Mean[(12/Pi)phaseds[[n]]],{n,1,id[1][[7]]}],id[1],"Mean phasse dev. (hr)",{-2.5,2.6}]

CalibRImagePlot[fftper,id[1],"Period (hr)",{21,25.1}]

CalibRImagePlot[Table[StandardDeviation[phaseds[[n]]],{n,1,id[1][[7]]}],id[1],"STD phasse dev. (hr)"]

border=TakeBorder[ptbl,id[1]];

core=TakeCore[ptbl,id[1]];

ShellPhasePlot[ptbl,phaseds,id[1]]

ShellAmplitudePlot[ptbl,phaseds,id[1]]

shell=TakeShell[ptbl,id[1]];

ftbl=ClusterLabel[ptbl][[All,2]];

env=GetAllEnvelope[sall,id[1]];

CLenv=Flatten[Table[env[[ptbl[[n]]]],{n,1,Length[ptbl]}],1];

RasterPlot[env,id[1],ptbl,CalibrationBarLabel->"Envelope",TimeLabel->"ExT",
MainLabel->"Ranked by mean luminescence",Clustered->False]

CLenv=Flatten[Table[Join[env[[border[[n]]]],env[[core[[n]]]]],{n,1,Length[ptbl]}],1];

RasterPlot[CLenv,id[1],ptbl,CalibrationBarLabel->"Envelope",TimeLabel->"ExT",
MainLabel->"Ranked by cluster",Clustered->True]

CLrfall=Flatten[Table[Join[rfall[[border[[n]]]],rfall[[core[[n]]]]],{n,1,Length[ptbl]}],1];

RasterPlot[CLrfall,id[1],ptbl,CalibrationBarLabel->"Normalized luminescence",TimeLabel->"ExT",
MainLabel->"Ranked by cluster (core->border)",Clustered->True]

RImageArray[rfall,id[1],MaxCycle->7,TimeLabel->"ExT",Color->1,Koma->2]

RImageArray[rlall,id[1],MaxCycle->7,TimeLabel->"ExT",Color->1,Koma->2]

RImageArray[rfall,id[1],MaxCycle->7,TimeLabel->"ExT",Color->1,Koma->3]

RImageArray[rlall,id[1],MaxCycle->7,TimeLabel->"ExT",Color->1,Koma->3]

extZero=Round[(24-id[1][[6]])]*4

IDID[6]

id[1][[1]]

ListPlot[Mean[Sin[abrac]][[Range[1,extZero+7]]]]

CalibRayleighPlot[phaseds,ptbl,0,0]

CalibRayleighPlot[abrac,ptbl,0,0]

Show[PlotClusterTimeSeries[sall,ptbl,id[1],TimeLabel->"ExT",Backdrop->True],
PlotRange->{{0,24*(60/id[1][[1]])+id[1][[2]]},{-Max[{Max[sall[[1]]],Abs[Min[sall[[1]]]]}],Max[{Max[sall[[1]]],Abs[Min[sall[[1]]]]}]}},AspectRatio->0.5]

Show[PlotClusterTimeSeries[rfall,ptbl,id[1],TimeLabel->"ExT",Backdrop->True],
PlotRange->{{0,24*(60/id[1][[1]])+id[1][[2]]},{-1.05,1.05}},AspectRatio->0.5]

PerPhaPlot[fftper,phaseds,ptbl,id[1]]

PeriodPhasePlot[fftper,pha,id[1]]

Show[PlotMeanTimeSeries[rlall,ptbl,id[1],TimeLabel->"ExT",Backdrop->True],
AspectRatio->0.5,PlotRange->{{0,24*(60/id[1][[1]])+id[1][[2]]},{-1.2,1.2}}]

Show[PlotMeanTimeSeries[rlall,ptbl,id[1],TimeLabel->"ZT",Backdrop->True],
AspectRatio->0.5,PlotRange->{{0,24*(60/id[1][[1]])+id[1][[2]]},{-1.2,1.2}}]

Show[PlotMeanTimeSeries[rlall,ptbl,id[1],Backdrop->False],
AspectRatio->0.5,PlotRange->{{0,144*2},{-1.2,1.2}}]

DoublePlot[rlall,ptbl,id[1]]