Package and Notebook code

.gitattributes

@@ -1,17 +0,0 @@
-# Auto detect text files and perform LF normalization
-* text=auto
-
-# Custom for Visual Studio
-*.cs     diff=csharp
-
-# Standard to msysgit
-*.doc	 diff=astextplain
-*.DOC	 diff=astextplain
-*.docx diff=astextplain
-*.DOCX diff=astextplain
-*.dot  diff=astextplain
-*.DOT  diff=astextplain
-*.pdf  diff=astextplain
-*.PDF	 diff=astextplain
-*.rtf	 diff=astextplain
-*.RTF	 diff=astextplain

.gitignore

@@ -1,47 +0,0 @@
-# Windows image file caches
-Thumbs.db
-ehthumbs.db
-
-# Folder config file
-Desktop.ini
-
-# Recycle Bin used on file shares
-$RECYCLE.BIN/
-
-# Windows Installer files
-*.cab
-*.msi
-*.msm
-*.msp
-
-# Windows shortcuts
-*.lnk
-
-# =========================
-# Operating System Files
-# =========================
-
-# OSX
-# =========================
-
-.DS_Store
-.AppleDouble
-.LSOverride
-
-# Thumbnails
-._*
-
-# Files that might appear in the root of a volume
-.DocumentRevisions-V100
-.fseventsd
-.Spotlight-V100
-.TemporaryItems
-.Trashes
-.VolumeIcon.icns
-
-# Directories potentially created on remote AFP share
-.AppleDB
-.AppleDesktop
-Network Trash Folder
-Temporary Items
-.apdisk

GenChaosGame.m

@@ -0,0 +1,49 @@
+(* ::Package:: *)
+
+BeginPackage["GenChaosGame`"]
+
+
+sparseFCGR::usage="sparseCGR[seq_,alphabet_List,k_] returns the FCGR plot of a genome as a sparse matrix. assumes linear, non-circular dna string"
+FCGR::usage="FCGR[seq_,alphabet_List,k_] returns the grayscale Chaos Game Representation plot of a DNA sequence in a unit square in 2D using 2^k resolution.assumes linear, non-circular dna string"
+
+Begin["`Private`"]
+
+sparseFCGR[seq_,alphabet_List,k_]:=
+Module[{chars,shifts,helpfunc1,helpfunc2,helpfunc3,pts,arraypts},
+chars=StringCases[seq,alphabet[[1]]|alphabet[[2]]|alphabet[[3]]|alphabet[[4]]];
+shifts=chars/.{alphabet[[1]]->{0.,0.},alphabet[[2]]->{0.,1.0*2^k},alphabet[[3]]->{1.0*2^k,1.0*2^k},alphabet[[4]]->{1.0*2^k,0.}};
+helpfunc1=Compile[{{a,_Real,1}},FoldList[IntegerPart[(#+#2)/2]&,2^(k-1),a]];
+helpfunc2=Compile[{{a,_Integer,1}},Map[2^k-#&,a]];
+helpfunc3=Compile[{{a,_Integer,1}},Map[1+#&,a]];
+pts=Transpose[helpfunc1/@Transpose[shifts]];
+pts=Drop[pts,k];
+arraypts=Transpose[{helpfunc2[Transpose[pts][[2]]],helpfunc3[Transpose[pts][[1]]]}];
+Return[{Length[pts],Tally[arraypts]}];
+];
+
+
+
+FCGR[seq_,alphabet_List,k_]:=
+Module[{resimg,chars,shifts,helpfunc1,helpfunc2,helpfunc3,pts,arraypts,setgrayscale},
+setgrayscale[{a_,b_}]:=Module[{oldval,newval},oldval=resimg[[a,b]];newval=oldval+1;resimg[[a,b]]=newval;];
+resimg=Table[0,{i,1,2^k},{j,1,2^k}];
+chars=StringCases[seq,alphabet[[1]]|alphabet[[2]]|alphabet[[3]]|alphabet[[4]]];
+shifts=chars/.{alphabet[[1]]->{0.,0.},alphabet[[2]]->{0.,1.0*2^k},alphabet[[3]]->{1.0*2^k,1.0*2^k},alphabet[[4]]->{1.0*2^k,0.}};
+helpfunc1=Compile[{{a,_Real,1}},FoldList[IntegerPart[(#+#2)/2]&,2^(k-1),a]];
+helpfunc2=Compile[{{a,_Integer,1}},Map[2^k-#&,a]];
+helpfunc3=Compile[{{a,_Integer,1}},Map[1+#&,a]];
+pts=Transpose[helpfunc1/@Transpose[shifts]];
+pts=Drop[pts,k];
+arraypts=Transpose[{helpfunc2[Transpose[pts][[2]]],helpfunc3[Transpose[pts][[1]]]}];
+Map[setgrayscale[{#[[1]],#[[2]]}]&,arraypts];
+Return[{Length[pts],resimg}];
+];
+
+
+End[]
+
+EndPackage[]
+
+
+
+

GenDistanceFunctions.m

@@ -0,0 +1,50 @@
+(* ::Package:: *)
+
+BeginPackage["GenDistanceFunctions`"]
+
+SSIMmax::usage="SSIMmax[img1_,img2_,max_] computes SSIM distance"
+mydescriptor::usage="mydescriptor[img_,winsizes_,steps_,histbins_] = Returns the descriptor of a  matrix/image, given the windowsizes, windowsteps and histogram bins"
+ApproxInfoDist::usage="Computes ApproxInfoDist"
+
+
+Begin["`Private`"]
+
+SSIMmax[img1_,img2_,max_]:=Module[{w, c1, c2, m1, m2, m1sq, m2sq, m1m2, sigma1sq, sigma2sq,sigma12, ssimmap, mssim},
+w = GaussianMatrix[{5, 1.5}]; c1 = (0.01*max)^2; c2 = (0.03*max)^2;
+m1 = ListCorrelate[w, img1]; m2 = ListCorrelate[w, img2];
+m1sq = m1*m1; m2sq = m2*m2; m1m2 = m1*m2;
+sigma1sq = ListCorrelate[w, img1*img1] - m1sq;
+sigma2sq = ListCorrelate[w, img2*img2] - m2sq;
+sigma12 = ListCorrelate[w, img1*img2] - m1m2; 
+ssimmap = ((c1 + 2*m1m2)*(c2 + 2*sigma12))/((c1 + m1sq + m2sq)*(c2 + sigma1sq + sigma2sq));
+mssim = Mean[Mean[ssimmap]];
+Return[mssim]];
+
+
+mydescriptor[img_,winsizes_,steps_,histbins_]:=
+Module[{dim=Length[img],allwinds={},winds,i,j,i1,nowsize,nowstep},
+If[Length@winsizes!=Length@steps,Print["Error in winsizes and steps lists."];Return[]];
+For[i1=1,i1<=Length[winsizes],i1++,
+nowsize=winsizes[[i1]];nowstep=steps[[i1]];
+winds=Table[Take[img,{i,nowsize+i-1},{j,nowsize+j-1}],{i,1,dim-nowsize+1,nowstep},{j,1,dim-nowsize+1,nowstep}];
+AppendTo[allwinds,Flatten[Table[(Length/@BinLists[Sort[Flatten[winds[[i,j]]]],{histbins}])/(nowsize^2),{i,1,Length[winds]},{j,1,Length[winds]}]]];
+];
+Return[Flatten[allwinds]]
+];
+
+ApproxInfoDist[i_, j_] := 
+Module[{x, y, xy, yx, cgr1, cgr2},
+cgr1 = i; cgr2 = j;
+x = Total[Unitize[cgr1], 2]; y = Total[Unitize[cgr2], 2];
+xy = Total[Unitize[cgr1 + cgr2], 2];
+yx = Total[Unitize[cgr1 + cgr2], 2];(*Print[{x,y,xy,yx}];*)
+Return[N[(xy - y + yx - x)/xy]];
+];
+
+End[]
+
+EndPackage[]
+
+
+
+

GenMDSTools.m

@@ -0,0 +1,67 @@
+(* ::Package:: *)
+
+BeginPackage["GenMDSTools`"]
+
+getUpperHalf::usage="getUpperHalf[matr_List] accepts as input an nxn matrix and returns the upper half of it (excluding the diagonal), flattened"
+mds::usage="mds[delta_List,dimensions_,dbg_] gets as input the distances delta_{ij}, and performs MDS"
+
+Begin["`Private`"]
+
+getUpperHalf[matr_List]:=
+Module[{i,res={}},
+For[i=1,i<=Length[matr]-1,i++,AppendTo[res,Take[matr[[i]],{i+1,Length[matr]}]]];Return[Flatten[res]]
+];
+
+mds[delta_List,dimensions_,dbg_,dim1_:1,dim2_:2,dim3_:3]:= (* DEBUG options= -1: PRINT NOTHING 0: all - 1: eigenvals+sol+MDSError *)
+Module[{n,deltasq,deltatotals,sumOfDelta,bMatr,eigensys,diagMatr,solpts,dMatr,mdserrors,pointstress,nofmdserrors,pos,step,finalpos},
+Monitor[
+n=Length[delta];
+deltasq=delta^2;
+step=1;
+deltatotals=Plus@@deltasq;
+sumOfDelta=Plus@@deltatotals;
+step=2;
+bMatr=Table[-(1/2)*(deltasq[[i, j]]-(1/n)*deltatotals[[i]]-(1/n)*deltatotals[[j]]+(1/(n^2))*sumOfDelta), {i, 1, n}, {j, 1, n}];
+(*jMatr=IdentityMatrix[Length[delta]]-(1/Length[delta])*Table[1,{k,1,Length[delta]},{l,1,Length[delta]}];*)
+(*bMatr=-1/2*jMatr.deltasq.jMatr;*)
+step=3;
+eigensys=N[Eigensystem[N[bMatr,20],Min[10,n]],20];  (*eigensys=N[Eigensystem[N[bMatr,15]],15];*)
+Print["10_First_Eigenvalues= ",eigensys[[1]]];
+step=4;
+pos=Select[Transpose[eigensys],#[[1]]>0&];
+If[dimensions==2,finalpos={pos[[dim1]],pos[[dim2]]};];
+If[dimensions==3,finalpos={pos[[dim1]],pos[[dim2]],pos[[dim3]]};];
+pos=Transpose[finalpos]; (* backward compatibility *)
+step=5;
+diagMatr=Sqrt[DiagonalMatrix[pos[[1]]]];  (* "Dimensions" largest positive eigenvalues *)
+step=6;
+solpts=Transpose[pos[[2]]].diagMatr;
+step=7;
+dMatr={};
+step=8;
+mdserrors=0;
+step=9;
+nofmdserrors=0;
+step=10;
+pointstress={};
+
+(* ALL dbg *)
+If[dbg==0,Print["delta=",delta//MatrixForm]];
+If[dbg==0,Print["bMatr=",N[bMatr,15]//MatrixForm]];
+If[dbg==0,Print["eigensys=",eigensys]];
+If[0<=dbg<=1,Print["eigenvals=",eigensys[[1]]]];
+If[dbg==0,Print["diagMatr=",diagMatr]];
+If[0<=dbg<=1,Print["sol=",solpts]];
+If[dbg==0,Print["dMatr=",dMatr//MatrixForm]];
+If[0<=dbg<=1,Print[mdserrors//MatrixForm]];
+If[dbg==0,Print["pointStress=",pointstress//MatrixForm]];
+If[0<=dbg<=1,Print[nofmdserrors//MatrixForm]];
+(* RETURN *)
+Return[{solpts,mdserrors,Select[eigensys[[1]],#>0&],pointstress,nofmdserrors}]
+,{step,{i,j},{i1,i2}}]];
+
+
+
+End[]
+
+EndPackage[]