tMake[table_]:=Module[{t,colheaders,rowheaders,colindex,rowindex,data},(*The (1,1) element of the input table is the row-index label*)t["row-index-type"]=table[[1,1]];
(*Extract the column and row headers*)colheaders=Rest[table[[1]]];
rowheaders=Rest[table[[All,1]]];
(*Place them in the tabular object*)t["column-headers"]=colheaders;
t["row-headers"]=rowheaders;
(*Compute the indexing functions*)For[i=1,i≤Length[colheaders],i=i+1,colindex[colheaders[[i]]]=i];
For[i=1,i≤Length[rowheaders],i=i+1,rowindex[rowheaders[[i]]]=i];
(*Install the indexing functions*)t["column-indexer"]=colindex;
t["row-indexer"]=rowindex;
(*Extract the actual data from the input matrix and install it*)data=table[[2;;-1,2;;-1]];
t["data"]=data;
(*Return the tabular object*)t];

tMake[t_,data_]:=Module[{full},(*Create a new table with the basic table constructor*)tMake[tFull[t,data]]];

tMake[rowIndexType_,rowHeaders_,columnHeaders_,data_]:=Module[{},tMake[Transpose[Prepend[Transpose[Prepend[data,columnHeaders]],Prepend[rowHeaders,rowIndexType]]]]];

tMake[rowIndexType_,rowHeaders_,columnHeaders_]:=Module[{zeroes},zeroes=Table[0,{i,rowHeaders},{j,columnHeaders}];
tMake[rowIndexType,rowHeaders,columnHeaders,zeroes]];

tColumn[t_,j_Integer]:=Module[{},t["data"][[All,j]]];

tColumn[t_,j_String]:=Module[{},t["data"][[All,t["column-indexer"][j]]]];

tRow[t_,i_Integer]:=Module[{},t["data"][[i]]];

tRow[t_,i_String]:=Module[{},t["data"][[All,t["row-indexer"][i]]]];

tRowIndexType[t_]:=Module[{},t["row-index-type"]];

tItem[t_,i_Integer,j_Integer]:=Module[{},t["data"][[i,j]]];

tSetItem[t_,i_Integer,j_Integer,v_]:=Module[{},t["data"]=ReplacePart[t["data"],{i,j}->v];];

tItem[t_,i_String,j_String]:=Module[{},t["data"][[t["row-indexer"][i],t["column-indexer"][j]]]];

tSetItem[t_,i_String,j_String,v_]:=Module[{},t["data"]=ReplacePart[t["data"],{t["row-indexer"][i],t["column-indexer"][j]}->v];];

tColumnHeaders[t_]:=Module[{},t["column-headers"]];

tRowHeaders[t_]:=Module[{},t["row-headers"]];

tData[t_]:=Module[{},t["data"]];

tFull[t_,m_]:=Module[{},Transpose[Prepend[Transpose[Prepend[m,t["column-headers"]]],Join[{t["row-index-type"]},t["row-headers"]]]]];

tFull[t_]:=Module[{},tFull[t,t["data"]]];

tTranspose[t_]:=Module[{},tMake[Transpose[tFull[t]]]];

tMapItems[t_,f_]:=Module[{newdata},(*Map the function to each data item*)newdata=Map[f,t["data"],{2}];
(*Construct a new table based on the new data*)tMake[t,newdata]];

tMapRows[t_,f_]:=Module[{newdata},(*Map the function to the rows of the data*)newdata=Map[f,t["data"]];
(*Construct a new table based on the new data*)tMake[t,newdata]];

tMapColumns[t_,f_]:=Module[{newdata},(*Map the function to the rows of the data*)newdata=Transpose[Map[f,Transpose[t["data"]]]];
(*Construct a new table based on the new data*)tMake[t,newdata]];

tMapMatrix[t_,f_]:=Module[{newdata},(*Call the function on the data matrix as a whole.*)newdata=f[t["data"]];
(*Construct a new table based on the new data*)tMake[t,newdata]];

tFilterRowData[t_,f_]:=Module[{rowdata},(*The row data is all but the top row of the full table representation*)rowdata=Rest[tFull[t]];
(*Filter the rows by applying the function to the data portion of each*)rowdata=Select[rowdata,f[Rest[#]]&];
(*Stick the correct column headers onto the row data,and recreate a table from it.*)PrependTo[rowdata,Prepend[t["column-headers"],t["row-index-type"]]];
tMake[rowdata]];

tFilterColumnData[t_,f_]:=Module[{coldata},(*The column data is all but the first column of the full table representation (taken here as the rows of the transpose of the full representation)*)coldata=Rest[Transpose[tFull[t]]];
(*Filter the columns by applying the predicate*)coldata=Select[coldata,f[Rest[#]]&];
(*Stick the row headers onto the column data,untranspose,and recreate a table from it.*)PrependTo[coldata,Prepend[t["row-headers"],t["row-index-type"]]];
tMake[Transpose[coldata]]];

tFilterRowHeaders[t_,f_]:=Module[{rowdata},(*This function repeats the actions of tFilterRowData[],but applies the predicate to the row header instead of the data*)rowdata=Rest[tFull[t]];
rowdata=Select[rowdata,f[First[#]]&];
PrependTo[rowdata,Prepend[t["column-headers"],t["row-index-type"]]];
tMake[rowdata]];

tFilterColumnHeaders[t_,f_]:=Module[{coldata},(*This function repeats the actions of tFilterColumnData[],but applies the predicate to the row header instead of the data*)coldata=Rest[Transpose[tFull[t]]];
coldata=Select[coldata,f[First[#]]&];
PrependTo[coldata,Prepend[t["row-headers"],t["row-index-type"]]];
tMake[Transpose[coldata]]];

(*Some colors*)muteYellow=RGBColor[0.527,0.527,0.0];
muteRed=RGBColor[0.847,0.057,0.057];
muteBlue=RGBColor[0.316,0.316,0.991];
muteGreen=RGBColor[0.0,0.592,0.0];
mutePurple=RGBColor[0.718,0.0,0.718];
muterPurple=RGBColor[0.574,0.112,0.686];

(*Cell cycle classes,and associated colors*)
classes={"None","M/G1","G1","S",    "S/G2","G2/M", "Other"};
colors={Black,     Yellow,Green,Blue,Red,        Orange, Purple};

classificationColor["None"] = {0,0,0};
classificationColor["M/G1"] = {1,1,0};
classificationColor["G1"] = {0,0.5,0};
classificationColor["S"] = {0,0,1};
classificationColor["S/G2"] = {1,0,0};
classificationColor["G2/M"] = {1, 0.5, 0};
classificationColor["Other"] = {0.5, 0, 1};

(* Function for creating a hatching pattern texture *)
texture[low_, high_] := {Table[(0.5+0.5*Sin[2*N[Pi]*20*(i+j)/256])*(high - low) + low,{i,256}, {j,256}]};

(* Function for creating a label for the final display *)
(*label[text_, width_, height_, texture_] :=
Module[{pts},
pts = {{0,0}, {width, 0}, {width, height}, {0, height}};
Graphics[{Texture[texture], (*EdgeForm[Thick],*) Polygon[pts, VertexTextureCoordinates -> pts],
		(*Inset[text, {width/2, height + 0.15}, BaseStyle -> {FontFamily -> "Courier"}]},*)
		Inset[text, Scaled[{0, 0.30}, {width/2, height}], BaseStyle -> {FontFamily -> "Courier"}]},
		Background -> GrayLevel[204/255]]
	];*)

label[text_, width_, height_, texture_] :=
Module[{pts},
pts = {{0,0}, {width, 0}, {width, height}, {0, height}};
Graphics[{Texture[texture], (*EdgeForm[Thick],*) Polygon[pts, VertexTextureCoordinates -> pts]}, ImageSize -> Automatic]
	];

horzLabel[text_, width_, height_, texture_] :=
Module[{pts},
pts = {{0,0}, {height, 0}, {height, width}, {0, width}};
Graphics[{Texture[texture], Polygon[pts, VertexTextureCoordinates -> pts]}, ImageSize -> Automatic]
];

(*Simple linear interpolation between two values.*)
linearInterpolate[a_,b_,t_]:=Module[{},t*(b-a)+a];

(*Interpolate colors,with desaturation at 0*)
colorSmooth[lowA_RGBColor,highA_RGBColor,valA_]:=Module[{low,high,val,col},(*Make the color args into lists*)low=Level[lowA,1];
high=Level[highA,1];
(*Transform the interpolation value from[-1,1] to[0,1]*)val=0.5*(valA+1.0);
(*Interpolate the color*)col=Table[linearInterpolate[low[[i]],high[[i]],val],{i,3}];
(*Fade the color toward white,depending on the closeness of the untransformed interpolation value to 0*)col=Table[linearInterpolate[col[[i]],1.0,1.0-Abs[val]],{i,3}];
(*Convert the vector into an RGBColor*)Apply[RGBColor,col]];

(*Use just two colors,with desaturation at 0*)
colorSharp[low_RGBColor,high_RGBColor,val_]:=Module[{col},(*Interpolate from white to one of the diverging colors,inverting the value if necessary*)col=If[val>0,Table[linearInterpolate[1.0,high[[i]],val],{i,3}],Table[linearInterpolate[1.0,low[[i]],-val],{i,3}]];
(*Turn the interpolated vector into a color*)Apply[RGBColor,col]];

(*Graphics primitive for ellipse drawing*)
ellipse[a_,b_,c_,theta_]:=Module[{npoints,x,y,ang,total,alpha,rad,xp,yp,p},npoints=100;
(*Construct a list of angles around the circle*)total=2*Pi;
ang=Table[x,{x,0,total,total/(npoints+1)}];
(*Compute the ellipse points*)(*Print["cos sin"];*)x=a Cos[ang];
y=b Sin[ang];
p=Transpose[{x,y}];
(*Print["alpha"];*)alpha={};
For[i=1,i≤Length[p],i++,AppendTo[alpha,vectorangle[p[[i]]]]];
(*Print["radius"];*)rad=Sqrt[x^2+y^2];
(*Print["cos sin 2"];*)xp=rad Cos[alpha+theta]+c[[1]];
yp=rad Sin[alpha+theta]+c[[2]];
(*Draw the ellipse*)p=Transpose[{xp,yp}];
Line[p]];

vectorangle[p_]:=Module[{a},a=ArcTan[p[[1]],Chop[p[[2]],10^-9]];
If[a<0,a=a+2Pi,Null];
If[Head[a]==Head[I],Print["NOOOO"];
Print[p];
Print[Chop[p]];
Print[a],Null];
a];

assemblePlot[plot_,ellipse_,correlation_,negColor_RGBColor,posColor_RGBColor]:=Module[{},(*Print[colorSharp[negColor,posColor,correlation]];*)Show[plot,Graphics[ellipse],AspectRatio->1,Background->colorSharp[negColor,posColor,correlation]]];

(*Given two input tables and a table describing the mappings to cell cycle classifications,creates a list of lists of points,as well as the corresponding cell cycle classifications (for use in a multicolored scatterplot*)
cellcycleClassPoints[expr_,prot_,cellcycle_]:=Module[{classBreakdown,exprBreakdown,protBreakdown,c},(*Filter the cellcycle list multiple times by class*)classBreakdown=Table[tFilterRowData[cellcycle,#[[1]]==c&],{c,classes}];
(*The input tables are now row filtered by membership of their row headers in each class.*)exprBreakdown=Table[tFilterRowHeaders[expr,MemberQ[tRowHeaders[c],#]&],{c,classBreakdown}];
protBreakdown=Table[tFilterRowHeaders[prot,MemberQ[tRowHeaders[c],#]&],{c,classBreakdown}];
(*Tag the tables with the labels and return them.*)Transpose[{classes,protBreakdown,exprBreakdown}]];

(*Computes the first two principal components of a set of points*)
pca[x_,y_]:=Module[{obs,ncols},(*Construct an observation matrix from two input vectors*)obs=Transpose[{x-Mean[x],y-Mean[y]}];
(*Compute the covariance matrix*)ncols=Length[x];
obs=Transpose[obs].obs/(ncols-1);
(*Compute the eigendecomposition of the covariance matrix-this contains the first two variances and modes of the PCA*)Eigensystem[obs]];

(*Annihilates the first singular mode of a matrix*)
svdCenterMatrix[m_]:=Module[{u,s,v},{u,s,v}=SingularValueDecomposition[m];
s[[1,1]]=0;
u.s.Transpose[v]];

(*Reconstructs matrix m from its n largest singular modes*)
subspaceReduceMatrix[m_,n_]:=Module[{u,s,v},{u,s,v}=SingularValueDecomposition[m,n];
u.s.Transpose[v]];

(*Projects the columns of one matrix onto the column space of another,in pseudoinverse fashion*)
pseudoinverseProjectionMatrix[basis_,data_,k_]:=Module[{u,s,v},(*SVD of the basis matrix*){u,s,v}=If[k==0,SingularValueDecomposition[basis],SingularValueDecomposition[basis,k]];
(*Pseudoinverse projection of the data matrix against the basis*)u.Transpose[u].data];

pseudoinverseProjectionMatrix[basis_,data_]:=pseudoinverseProjectionMatrix[basis,data,0];

(**************************************)

(*Read an xls file in,and construct a table from it.*)
readXLS[file_,i_]:=Module[{},tMake[Import[file][[i]]]];

readXLS[file_]:=readXLS[file,1];

readTable[file_]:=Module[{},tMake[Import[file,"TSV"]]];

(*These functions eliminate rows/columns of a table containing string data.*)
eliminateNonnumericColumns[t_]:=Module[{n},n=Length[tRowHeaders[t]];
tFilterColumnData[t,Length[Select[#,Function[{x},Head[x]===String]]]==0&]];

(*The following functions are used to select rows/columns of a pvalue file based on how many of the values fall above a threshold*)
valueCountFilterColumns[pval_,threshold_,frac_]:=Module[{n,passedItems},(*Calculate the minimum number of pvalues we need to see above the threshold in order to keep a protein*)n=Ceiling[Length[tColumnHeaders[pval]]*frac];
(*The function argument says "the number of items left
              after only keeping those items greater than threshold is
              greater than n"*)tFilterColumnData[pval,Length[Select[#,Function[x,x>threshold]]]>n]];

valueCountFilterRows[pval_,threshold_,frac_]:=Module[{},valueCountFilterColumns[tTranspose[pval],threshold,frac]];

(*Take "raw" basis and data files and transform them into a form suitable for correlation-seeking*)
process[basisArg_,dataArg_,pipMode_String,pipColumnsArg_Integer]:=Module[{basis,data,pipColumns,dataPIP},basis=basisArg;
data=dataArg;
(*If the pipColumns argument is 0,this actually means to use all of the columns during PIP and subspace reduction*)pipColumns=If[pipColumnsArg==0,Dimensions[tData[basis]][[2]],pipColumnsArg];
(*Remove columns containing non-numeric data from both datasets.*)Print["eliminate"];
basis=eliminateNonnumericColumns[basis];
data=eliminateNonnumericColumns[data];
(*Compute the row intersection of the two tables.*)Print["intersect"];
{basis,data}=rowIntersect[basis,data];
(*SVD center the data sets*)Print["center"];
basis=svdCenter[basis];
data=svdCenter[data];
(*Perform PIP of data against basis*)Print["PIP"];
dataPIP=pseudoinverseProjection[basis,data,pipColumns];
(*Take other actions for different PIP modes*)Print["PIP-followon"];
Switch[pipMode,(*Do nothing in the "regular" case*)"regular",Null,(*Perform the basis PIP as well (against the original data)*)"mutual",basis=pseudoinverseProjection[data,basis,pipColumns],(*Perform subspace reduction on the basis matrix*)"subspace",basis=subspaceReduce[basis,pipColumns],_,"Error"];
(*Let "data" become its PIP*)data=dataPIP;
(*Center the columns of both matrices*)Print["column-center"];
basis=columnCenter[basis];
data=columnCenter[data];
(*Normalize the columns of both matrices*)Print["column-normalize"];
basis=columnNormalize[basis];
data=columnNormalize[data];
(*Return the transformed basis and data*){basis,data}];

(*Compute the correlations between protein bindings and gene expression levels*)
correlation[basis_,data_]:=correlation[basis,data,Correlation];

correlation[basis_,data_,corrFunc_]:=Module[{corrdata},(*Compute correlations of each protein and each time point.*)corrdata=Table[N[corrFunc[tColumn[data,p],tColumn[basis,t]]],{p,tColumnHeaders[data]},{t,tColumnHeaders[basis]}];
(*Construct a table of the correlation values*)tMake["",tColumnHeaders[data],tColumnHeaders[basis],corrdata]];

regressionCorrelation[regr_]:=Module[{corr,sd,topmode},corr=Table[sd=tItem[regr,p,t][[1]];
topmode=tItem[regr,p,t][[2]];
Sign[topmode[[1,1]]]*Sign[topmode[[1,2]]]*(1.0-sd[[2]]/sd[[1]]),{p,tRowHeaders[regr]},{t,tColumnHeaders[regr]}];
tMake["",tRowHeaders[regr],tColumnHeaders[regr],corr]];

(*Compute the regression lines using PCA between the protein bindings and gene expression levels*)
regression[basis_,data_]:=Module[{reg},(*Compute regressions for each protein and time point*)reg=Table[pca[tColumn[data,p],tColumn[basis,t]],{p,tColumnHeaders[data]},{t,tColumnHeaders[basis]}];
(*Construct a table of the regression values*)tMake["",tColumnHeaders[data],tColumnHeaders[basis],reg]];

plot[basis_,data_,regressions_,range_,ticks_]:=
Module[{plots,x,y,angle,ellipses,frame, size},
(*frame=If[ticks,{{{range[[1,1]],0,range[[1,2]]},None},{{range[[2,1]],0,range[[2,2]]},None}},None];*)
frame=If[ticks,{{range[[1,1]]*0.67,0,range[[1,2]]*0.67},{range[[2,1]]*0.67,0,range[[2,2]]*0.67}},None];

size = If[ticks, 0.2, 1];
(*Form data for a table of scatter plots*)plots=Table[(*Get the points to plot*)x=tColumn[data,p];
y=tColumn[basis,t];
ListPlot[Transpose[{x,y}],(*PlotRange->{{-0.1,0.1},{-0.1,0.1}},*)PlotRange->range,(*PlotRangePadding -> None,*) PlotStyle->GrayLevel[0.2],AspectRatio->1.0,PlotMarkers->{Automatic, 4},Frame->False,(*FrameTicks->{{Automatic,None},{Automatic,None}}*)Ticks->frame, ImageSize -> Automatic
(* ,Ticks->{{-0.1,0.1},{-0.1,0.1}}*)],{p,tColumnHeaders[data]},{t,tColumnHeaders[basis]}];
(*Form data for a table of ellipses*)ellipses=Table[(*Get the PCA modes in order to create an ellipse*){sd,mode}=tItem[regressions,p,t];
(*Compute an ellipse (along with the angle by which to rotate it)*)angle=vectorangle[mode[[1]]];
(*Print["ellipse "<>p<>t];*)ellipse[3*Sqrt[sd[[1]]],3*Sqrt[sd[[2]]],{Mean[x],Mean[y]},angle],{p,tColumnHeaders[data]},{t,tColumnHeaders[basis]}];
{tMake["",tColumnHeaders[data],tColumnHeaders[basis],plots],tMake["",tColumnHeaders[data],tColumnHeaders[basis],ellipses]}];

(*coloronlyPlot[correlations_,negColor_,posColor_]:=
Module[{plots},
plots=Table[Graphics[{colorSharp[negColor,posColor,tItem[correlations,p,t]],Rectangle[], ImageSize -> Automatic}],
{p,tRowHeaders[correlations]},{t,tColumnHeaders[correlations]}];
tMake["",tRowHeaders[correlations],tColumnHeaders[correlations],plots]];
*)

coloronlyPlot[correlations_,negColor_,posColor_]:=
Module[{plots},
plots=Table[Show[ListPlot[{{1,1}}, PlotRange -> {{0,0},{0.1,0.1}}, AspectRatio -> 1.0, Frame -> False, Axes -> False, ImageSize -> Automatic], Background -> colorSharp[negColor, posColor, tItem[correlations,p,t]]],
{p,tRowHeaders[correlations]},{t,tColumnHeaders[correlations]}];
tMake["",tRowHeaders[correlations],tColumnHeaders[correlations],plots]];

colorPlot[breakdown_,prots_,times_,range_]:=Module[{points,plots,prot,time},(*Collect several lists of points from the breakdown.*)points=Table[prot=breakdown[[i]][[2]];
time=breakdown[[i]][[3]];
Transpose[{tColumn[prot,p],tColumn[time,t]}],{p,prots},{t,times},{i,Length[breakdown]}];
(*Construct the plots*)plots=Table[ListPlot[points[[i,j]],PlotRange->range,PlotStyle->Rest[colors],AspectRatio->1.0,PlotMarkers->None,Frame->True,(*FrameTicks->{{Automatic,None},{Automatic,None}}],*)FrameTicks->None],{i,Dimensions[points][[1]]},{j,Dimensions[points][[2]]}];
(*Make a table of the plots and return it.*){points,tMake["",prots,times,plots]}];

assemble[plots_,ellipses_,correlations_,negColor_,posColor_]:=Module[{data},data=Table[assemblePlot[tItem[plots,p,t],tItem[ellipses,p,t],tItem[correlations,p,t],negColor,posColor],{p,tRowHeaders[plots]},{t,tColumnHeaders[plots]}];
tMake["",tRowHeaders[plots],tColumnHeaders[plots],data]];

export[plots_]:=Module[{},Table[Export[p<>"_"<>t<>".png",tItem[plots,p,t]],{p,tRowHeaders[plots]},{t,tColumnHeaders[plots]}];
Null];

(**************************************)

(*Compute the row-wise intersection of two tables-that is,throw out all rows with headers that do not occur in both tables*)
rowIntersect[t1_,t2_]:=Module[{rows1,rows2,common,i1,i2},(*Grab the row headers from both tables,and get their intersection*)rows1=tRowHeaders[t1];
rows2=tRowHeaders[t2];
common=Intersection[rows1,rows2];
(*Filter the row headers of each table by membership in the intersection*)i1=tFilterRowHeaders[t1,MemberQ[common,#]&];
i2=tFilterRowHeaders[t2,MemberQ[common,#]&];
{i1,i2}];

(*Convert the input table to "ratio format"-that is,each row's entries are divided by that row's mean.*)
ratioFormat[t_]:=Module[{},tMapRows[t,#/Mean[#]&]];

(*Center the rows of a table*)
rowCenter[t_]:=Module[{},tMapRows[t,#-Mean[#]&]];

(*Get the means of the columns of a table*)
columnMeans[t_]:=Module[{},Table[Mean[tColumn[t,i]],{i,Length[tColumnHeaders[t]]}]];

(*Center the columns of a table*)
columnCenter[t_]:=Module[{},tMapColumns[t,#-Mean[#]&]];

(*Perform "SVD centering" on a table*)
svdCenter[t_]:=Module[{},tMapMatrix[t,svdCenterMatrix]];

(*Subspace reduction on a table*)
subspaceReduce[t_,n_]:=Module[{},tMapMatrix[t,subspaceReduceMatrix[#,n]&]];

(*Pseudoinverse projection on tables*)
pseudoinverseProjection[basis_,data_,k_]:=Module[{},tMapMatrix[data,pseudoinverseProjectionMatrix[tData[basis],#,k]&]];

pseudoinverseProjection[basis_,data_]:=Module[{},tMapMatrix[data,pseudoinverseProjectionMatrix[tData[basis],#]&]];

(*Normalize the columns of a table*)
columnNormalize[t_]:=Module[{},tMapColumns[t,If[Chop[Norm[#]]==0,#,#/Norm[#]]&]];

classify[aArg_]:=Module[{a},
a = aArg;
While[a > 2*Pi, a = a - 2*Pi];

If[a≤Pi/2,Return["M/G1"],Null];
If[a≤Pi,Return["G1"],Null];
If[a≤Pi+(2/3)*(Pi/2),Return["S"],Null];
If[a≤3*Pi/2+(1/3)*(Pi/2),Return["S/G2"],Null];
If[a≤2*Pi,Return["G2/M"],Null];
Print[a];
Null];
SetAttributes[classify,Listable];

classifyS[a_] :=
Module[{g1,sg2, g2m, mg1, s, other},
g1 = {"MBP1", "SWI4", "SWI6", "14_min", "21_min", "28_min", "77_min", "84_min", "91_min"};
sg2 = {"FKH1", "42_min", "98_min"};
g2m = {"FKH2", "NDD1", "MCM1", "0_min", "49_min", "56_min", "63_min", "105_min", "112_min", "119_min"};
mg1 = {"ACE2", "SWI5", "7_min", "70_min"};
s = {"35_min"};
other = {"MCM3", "MCM4", "MCM7", "ORC1"};

If[MemberQ[g1, a], Return["G1"], Null];
If[MemberQ[sg2, a], Return["S/G2"], Null];
If[MemberQ[g2m, a], Return["G2/M"], Null];
If[MemberQ[mg1, a], Return["M/G1"], Null];
If[MemberQ[s, a], Return["S"], Null];
If[MemberQ[other, a], Return["Other"], Null];

Print["Error: "<>a<> " was not found in the known labels"]
];
SetAttributes[classifyS, Listable]

syntheticClassification[numGenes_]:=Module[{c},c=classify[Table[(a/numGenes)*2*Pi,{a,1,numGenes}]];
tMake["",Table["g"<>ToString[i],{i,1,numGenes}],{"synthetic_classification"},Transpose[{c}]]];

syntheticExpression[numGenes_,numTimes_]:=Module[{},Table[Sin[(a/numGenes-2*b/numTimes)*2*Pi],{a,1,numGenes},{b,1,numTimes}]/Sqrt[numGenes/2]];

syntheticExpressionTable[numGenes_,numProts_]:=Module[{colheaders,rowheaders},colheaders=Table["θ = "<>ToString[PaddedForm[N[p*2/numProts], {3,2}]]<>"π",{p,1,numProts}];
rowheaders=Table["g"<>ToString[g],{g,1,numGenes}];
data=syntheticExpression[numGenes,numProts];
tMake["",rowheaders,colheaders,N[data]]];

syntheticBinding[numGenes_,numProts_]:=Module[{},Table[Sin[(a/numGenes-b/numProts)*2*Pi],{a,1,numGenes},{b,1,numProts}]/Sqrt[numGenes / 2]];

syntheticBindingTable[numGenes_,numTimes_]:=Module[{colheaders,rowheaders},colheaders=Table["θ = "<>ToString[PaddedForm[N[p*2/numTimes] ,{3,2}] ]<>"π",{p,1,numTimes}];
rowheaders=Table["g"<>ToString[g],{g,1,numGenes}];
data=syntheticBinding[numGenes,numTimes];
tMake["",rowheaders,colheaders,N[data]]];

syntheticBindingClassification[numProts_] :=
	Module[{t,s,c},
		t = Table[(b/numProts)*2*Pi, {b, 1, numProts}];
		(*s = Table["("<>ToString[(b/numProts)*2, InputForm]<>")π", {b, 1, numProts}];*)
		s = Table[ToString[b*2]<>"π/"<>ToString[numProts], {b, 1, numProts}];
		c = classify[t];
		tMake["", s, {"synthetic_classification"}, Transpose[{c}]]   
	];

syntheticExpressionClassification[numTimes_] :=
	Module[{t,s,c},
		t = Table[(2*b/numTimes)*2*Pi, {b, 1, numTimes}];
		(*s = Table["("<>ToString[(2*b/numTimes)*2, InputForm]<>")π", {b, 1, numTimes}];*)
		s = Table[ToString[4*b]<>"π/"<>ToString[numTimes], {b, 1, numTimes}];
		c = classify[t];
		tMake["",s, {"synthetic_classification"}, Transpose[{c}]]   
	];

oneline[text_] := StringReplace[ToString[text,InputForm], {"Pi" -> "π", "*" -> "", "(" -> "", ")" -> ""}];

labels[classification_, size_] := 
Module[{t},
	Table[t = texture[{204/255,204/255,204/255}, classificationColor[tColumn[classification, 1][[i]] ] ];
		label[oneline[tRowHeaders[classification][[i]] ], size[[1]], size[[2]], t],
		{i, Length[tRowHeaders[classification]]}]
];

horzLabels[classification_, size_] :=
Module[{t},
	Table[t = texture[{204/255,204/255,204/255}, classificationColor[tColumn[classification, 1][[i]] ] ];
		horzLabel[oneline[tRowHeaders[classification][[i]] ], size[[1]], size[[2]], t],
		{i, Length[tRowHeaders[classification]]}]
];

(*assembleGrid[plots_, toplabels_, sidelabels_] := Transpose[Join[{Prepend[Table[Rotate[l, Pi/2], {l,sidelabels}], ""]},Transpose[Join[{toplabels}, tData[plots]]]]];*)
assembleGrid[plots_, toptext_, toplabels_, sidetext_, sidelabels_] :=
Transpose[
Join[
{Join[{""}, {""}, sidetext]},
{Join[{""}, {""}, sidelabels ]},
Transpose[Join[{toptext}, {toplabels}, tData[plots]]]
]
];

(* Set number of genes, proteins, and time points for creating synthetic data *)
numGenes = 20;
numProteins = 3;
numTimes = 7;
plotrange = {{-0.25, 0.25}, {-0.25, 0.25}};

(* Generate the synthetic data - flowchart phase "data processing" - no special processing needed due to design of synthetic data functions *)
binding = syntheticBindingTable[numGenes, numProteins];
expression = syntheticExpressionTable[numGenes, numTimes];
classification = syntheticClassification[numGenes];
bindingClassification = syntheticBindingClassification[numProteins];
expressionClassification = syntheticExpressionClassification[numTimes];

tFull[bindingClassification]

(* This calculation shows that the columns of the binding and expression matrices are normalized by design. *)
Table[Norm[tColumn[binding, i]], {i, Length[tColumnHeaders[binding]]}]
Table[Norm[tColumn[expression, i]], {i, Length[tColumnHeaders[expression]]}]

(* Display the data in grayscale raster images *)
scaledTable[t_] := Module[{tab},
						tab = tData[t];
						(tab - Min[tab]) / (Max[tab] - Min[tab])];
Grid[{{"Binding", "Expression"},{Graphics[Raster[scaledTable[binding]], AspectRatio -> GoldenRatio],Graphics[Raster[scaledTable[expression]], AspectRatio -> GoldenRatio]}}]

(* Compute the regressions using PCA - flowchart phase "Regression", step 1 *)
regr = regression[expression, binding];

(* Compute the statistical (Pearson) correlation - flowchart phase "Regression", step 2 *)
corr = correlation[expression, binding];
Grid[Chop[tFull[corr]], ItemSize -> Full]

(* This is the correlation function exactly as it appears in Jack's paper - note that it matches exactly with the Correlation[] function in mathematica *)
jack[x_, y_] :=
	Module[{sxy, sx, sy},
	sxy = Total[(x - Mean[x])*(y - Mean[y])] / (Length[x] - 1);
	sx = Sqrt[Total[(x - Mean[x])^2]/(Length[x] - 1)];
	sy = Sqrt[Total[(y - Mean[y])^2]/(Length[y] - 1)];
	sxy / (sx*sy)
];
pearsoncorr = correlation[expression, binding, jack];
Grid[Chop[tData[pearsoncorr] - tData[corr]]]

(* Compute data ellipses and scatterplots - flowchart phase "Plotting", steps 2 and 3 *)
{fullplots, ellipses} = plot[expression, binding, regr, plotrange, True];
{thumbplots, ellipses} = plot[expression, binding, regr, plotrange, False];

(* Create some labels for use in the plots *)
bindingLabels = horzLabels[bindingClassification, {1,0.25}];
bindingText = tRowHeaders[bindingClassification];
expressionLabels = labels[expressionClassification, {1, 0.25}];
expressionText = tRowHeaders[expressionClassification];

Grid[{{bindingLabels[[1]], expressionLabels[[1]]}}, ItemSize -> {{3,Automatic}, {Automatic, 3}}, Background -> Pink, Dividers -> All];

(* Draw the scatterplots overlaid with the ellipse, and colormapped to indicate the correlation value - flowchart phase "Plotting", step 1 *)
thumbassembledPurpleGreen = assemble[thumbplots, ellipses, corr, muterPurple, muteGreen];
(*thumbassembledRedBlue = assemble[thumbplots, ellipses, corr, muteRed, muteBlue];
thumbassembledRedGreen = assemble[thumbplots, ellipses, corr, muteRed, muteGreen];
thumbassembledYellowBlue = assemble[thumbplots, ellipses, corr, Yellow, muteBlue];
thumbassembled = {thumbassembledPurpleGreen, thumbassembledRedBlue, thumbassembledRedGreen, thumbassembledYellowBlue};*)

assembledPurpleGreen = assemble[fullplots, ellipses, corr, muterPurple, muteGreen];
(*assembledRedBlue = assemble[fullplots, ellipses, corr, muteRed, muteBlue];
assembledRedGreen = assemble[fullplots, ellipses, corr, muteRed, muteGreen];
assembledYellowBlue = assemble[fullplots, ellipses, corr, Yellow, muteBlue];
assembled = {assembledPurpleGreen, assembledRedBlue, assembledRedGreen, assembledYellowBlue};*)

coloronlyPurpleGreen = coloronlyPlot[corr, muterPurple, muteGreen];
(*coloronlyRedBlue = coloronlyPlot[corr, muteRed, muteBlue];
coloronlyRedGreen = coloronlyPlot[corr, muteRed, muteGreen];
coloronlyYellowBlue = coloronlyPlot[corr, Yellow, muteBlue];
coloronly = {coloronlyPurpleGreen, coloronlyRedBlue, coloronlyRedGreen, coloronlyYellowBlue};*)

(*allassembled = {thumbassembled, assembled, coloronly};*)

(* Display the assembled plots *)
(*DynamicModule[{colorscheme = 1, whichplots = 3},
Column[{Row[{RadioButtonBar[Dynamic[whichplots], {1->"frame labels off", 2->"frame labels on", 3->"color only"}]}],
	   Row[{RadioButtonBar[Dynamic[colorscheme], {1 -> "purple/green", 2 -> "red/blue", 3 -> "red/green", 4 -> "yellow/blue"}]}],
	   (*Dynamic[Grid[tFull[allassembled[[whichplots, colorscheme]] ],Spacings -> {0.1,0.1}, ItemSize -> Full*)

Dynamic[Panel[Grid[assembleGrid[allassembled[[whichplots, colorscheme]], expressionLabels, bindingLabels],(*Background -> GrayLevel[204/255],*) Spacings -> {0.1,0.1}, ItemSize -> Full  (*Scaled[1/Length[tColumnHeaders[allassembled[[whichplots, colorscheme]] ]]] *)], Background -> GrayLevel[204/255]] ]}]
];*)

makeGrid[plotsets_, whichColor_, whichPlot_, scale_, exprText_, exprLabels_, bindingText_, bindingLabels_, single_, i_, j_] :=
Module[{},
If[single,
Grid[assembleGrid[tMake["", {""}, {""}, {{tItem[plotsets[[whichColor, whichPlot]], i, j]}}], {exprText[[j]]}, {exprLabels[[j]]}, {bindingText[[i]]}, {bindingLabels[[i]]}],
Spacings -> {0.4, 0.4},
ItemSize -> {{Automatic, 3.7/scale+ 0.1, {15/scale}}}],
Grid[assembleGrid[plotsets[[whichColor, whichPlot]], exprText, exprLabels, bindingText, bindingLabels],
Spacings -> {0.4, 0.4},
ItemSize -> {{Automatic, 3.7/1+ 0.1, {15/1}}}]
]
];

Panel[makeGrid[{{coloronlyPurpleGreen, thumbassembledPurpleGreen, assembledPurpleGreen}}, 1, 1, 1.0, expressionText, expressionLabels, Table[Rotate[t,90Degree], {t,bindingText}], bindingLabels, False, 1, 1],
	Background -> GrayLevel[204/255]
];

(*DynamicModule[{colorscheme = 1, whichplots = 1, ij = {1,1}},
Column[
{Row[{RadioButtonBar[Dynamic[whichplots], {1->"color only", 2->"thumbnails", 3->"full plots"}], Slider2D[Dynamic[ij], {{1,1}, {numTimes, numProteins}, {1,1}}]}],
Row[{RadioButtonBar[Dynamic[colorscheme], {1->"purple/green"}]}],
Dynamic[Panel[makeGrid[{{coloronlyPurpleGreen, thumbassembledPurpleGreen, assembledPurpleGreen}}, colorscheme, whichplots, expressionText, expressionLabels, Table[Rotate[t,90Degree], {t,bindingText}], bindingLabels, whichplots == 3, ij[[2]], ij[[1]]],
	Background -> GrayLevel[204/255]
]]}
]
]*)

Manipulate[Panel[makeGrid[{{coloronlyPurpleGreen, thumbassembledPurpleGreen, assembledPurpleGreen}}, colorscheme, whichplots, 1.0, expressionText, expressionLabels, Table[Rotate[t,90Degree], {t,bindingText}], bindingLabels, whichplots == 3, i, j],
	Background -> GrayLevel[204/255]
],
{{whichplots, 1, "plot style"}, {1->"color only", 2->"thumbnails", 3->"full plots"}, ControlType -> RadioButtonBar},
{{colorscheme, 1, "color scheme"}, {1 -> "purple/green"}, ControlType -> RadioButtonBar},
{{i, 1, "zoom row"}, Table[x, {x, numProteins}], ControlType -> SetterBar},
{{j, 1, "zoom column"}, Table[x, {x, numTimes}], ControlType -> SetterBar},
SaveDefinitions -> True
]

p=Panel[Grid[assembleGrid[coloronlyPurpleGreen, expressionText, expressionLabels, Table[Rotate[t,90Degree], {t,bindingText}], bindingLabels],
		Spacings -> {0.4, 0.4},
	Dividers->None,
(*ItemSize -> {{Automatic, Scaled[0.26/Length[tColumnHeaders[coloronlyPurpleGreen]]], {Scaled[1/Length[tColumnHeaders[coloronlyPurpleGreen]]]}}}*)
(*ItemSize -> {{Automatic, 4.1, {16.4}}}*)
ItemSize -> {{Automatic, 3.7/1+ 0.1, {15/1}}}
(*ItemSize -> Full*)
	],
	Background -> GrayLevel[204/255]];

(*Export["test.pdf", p, ImageSize -> 750]*)

Panel[Grid[assembleGrid[thumbassembledPurpleGreen, expressionText, expressionLabels, Table[Rotate[t, 90Degree], {t, bindingText}], bindingLabels],
		Dividers -> None,
		Spacings -> {0.4, 0.4},
		(*ItemSize -> {{Automatic, Scaled[0.245/Length[tColumnHeaders[thumbassembledPurpleGreen]]], {Scaled[1/Length[tColumnHeaders[thumbassembledPurpleGreen]]]}}}*)
		ItemSize -> { {Automatic, 3.7/1+0.1, {15/1}}}
		(*ItemSize -> Automatic*)
	],
	Background -> GrayLevel[204/255]];

Panel[Grid[assembleGrid[assembledPurpleGreen, expressionText, expressionLabels, Table[Rotate[t, 90Degree], {t, bindingText}], bindingLabels],
		Dividers -> None,
		Spacings -> {0.4, 0.4},
		ItemSize -> { {Automatic, 3.7/1+0.1, {15/1}}}
	],
	Background -> GrayLevel[204/255]];

Panel[Grid[assembleGrid[assembledPurpleGreen, expressionText, expressionLabels, bindingText, bindingLabels],
		Spacings -> {0.1, 0.1},
		(*ItemSize ->{{Automatic, 3.5,{Full}}}*)
		ItemSize -> {{Automatic, Scaled[0.225*2/Length[tColumnHeaders[assembledPurpleGreen]]], {Scaled[2/Length[tColumnHeaders[assembledPurpleGreen]]]}}}
	],
	Background -> GrayLevel[204/255]];

(* Divide the points in each plot by which cell cycle they come from, and then create colored dot plots reflecting these cycles. *)
(*breakdown = Rest[cellcycleClassPoints[expression, binding, classification]];
{points, colorplots} = colorPlot[breakdown, tColumnHeaders[binding], tColumnHeaders[expression], plotrange];
assembledColor = assemble[colorplots, ellipses, corr, RGBColor[1,1,1], RGBColor[1,1,1]];*)

(*Grid[tFull[assembledColor], Spacings -> {0.1, 0.1}, ItemSize ->  Full  (*Scaled[1/Length[tColumnHeaders[assembledColor]]]*) ];*)

(* Create a unified image of the plots *)
(*grid = assembleGrid[assembledColor, expressionLabels, bindingLabels];
g = Grid[grid, Background -> GrayLevel[0.5], ItemSize -> Full];
*)

(* Export the grid as an image *)
(*Export["grid.png", g]*)

(* Repeat the process for real data *)
(* Read in the data *)
prot = readTable["/home/roni/work/genomics/data/original/Protein_Bind.txt"];
expr = readTable["/home/roni/work/genomics/data/original/Alpha_Factor.txt"];
protP = eliminateNonnumericColumns[prot];
exprP = eliminateNonnumericColumns[expr];

(* Process the data *)
{exprP, protP} = process[exprP, protP, "subspace", 4];

(* Compute the correlation *)
datacorr = correlation[exprP, protP];

(* Compute regressions *)
dataregr = regression[exprP, protP];

(* Compute the plots *)
{datafullplots, dataellipses} = plot[exprP, protP, dataregr, {{-0.1,0.1},{-0.1,0.1}}, True];
{datathumbplots, dataellipses} = plot[exprP, protP, dataregr, {{-0.1,0.1},{-0.1,0.1}}, False];

(* Assemble the plots *)
datathumbassembledPurpleGreen = assemble[datathumbplots, dataellipses, datacorr, mutePurple, muteGreen];
dataassembledPurpleGreen = assemble[datafullplots, dataellipses, datacorr, mutePurple, muteGreen];
datacoloronlyPurpleGreen = coloronlyPlot[datacorr, mutePurple, muteGreen];

(* Compute the labels *)
databindingText = tColumnHeaders[protP];
databindingClassification = tMake["", databindingText, {"class"}, Transpose[{classifyS[databindingText]}]];
databindingLabels = horzLabels[databindingClassification, {1,0.25}];

dataexpressionText = tColumnHeaders[exprP];
dataexpressionClassification = tMake["", dataexpressionText, {"class"}, Transpose[{classifyS[dataexpressionText]}]];
dataexpressionLabels = labels[dataexpressionClassification, {1, 0.25}];

tFull[databindingClassification]

Panel[Grid[assembleGrid[datacoloronlyPurpleGreen, dataexpressionText, dataexpressionLabels, databindingText, databindingLabels],
		Spacings -> {0.1, 0.1},
		ItemSize -> {{Automatic, Scaled[0.26/Length[tColumnHeaders[datacoloronlyPurpleGreen]]], {Scaled[1/Length[tColumnHeaders[datacoloronlyPurpleGreen]]]}}}
	],
	Background -> GrayLevel[204/255]];

Panel[Grid[assembleGrid[datathumbassembledPurpleGreen, dataexpressionText, dataexpressionLabels, databindingText, databindingLabels],
		Spacings -> {0.1, 0.1},
		ItemSize -> {{Automatic, Scaled[0.26/Length[tColumnHeaders[datathumbassembledPurpleGreen]]], {Scaled[1/Length[tColumnHeaders[datathumbassembledPurpleGreen]]]}}}
	],
	Background -> GrayLevel[204/255]];

Panel[Grid[assembleGrid[dataassembledPurpleGreen, dataexpressionText, dataexpressionLabels, databindingText, databindingLabels],
		Spacings -> {0.1, 0.1},
		ItemSize -> {{Automatic, Scaled[0.24*2/Length[tColumnHeaders[datathumbassembledPurpleGreen]]], {Scaled[2/Length[tColumnHeaders[datathumbassembledPurpleGreen]]]}}}
	],
	Background -> GrayLevel[204/255]];