Print["Updated: Feb. 6, 2001"]; 
(* Author: A.I. McLeod *)
(* Copyright (c), A.I. McLeod, 2001 *)
(* MViz *)
(* Summary:
Multivariate Visualization
*)



BeginPackage["MVIZ`"]
(*------------------------------------------------------*)
(*load standard packages*)

Needs["Statistics`NormalDistribution`"];
Needs["Graphics`Graphics`"];

(*----------------------------------------------------------------------*)
(*utilities*)

NumericVectorQ::usage="NumericVectorQ[x] returns true if x is a vector
and its elements are floating point numbers.";

NumberVectorQ::usage="NumberVectorQ[x] returns true if x is a vector
and its elements numbers (ie., rational, integer, floats).";

NumberListQ::usage="NumberListQ[x] is true if x is a list containing numbers.";

PowerTransformation::usage = 
"PowerTransformation[z, lambda] data z are transformed according to \
lambda. If lambda is 1, no transformation. If lambda is zero, logs are taken. \ 
Otherwise, z^lambda.";

RescaleDataMatrix::usage= 
"RescaleDataMatrix[x] each column is scaled independently to the 
range [0,1] after it has been centered -- mean corrected and divided
by its standard deviation.";
 
StandardizeDataMatrix::usage= 
"StandardizeDataMatrix[x] each column is rescaled to [0,1].";
 
StandardizeDataVector::usage=
"StandardizeDataVector[x] the vector is rescaled to [0,1].";
 
CenterDataMatrix::usage= 
"CenterDataMatrix[x] each column is rescaled to [0,1].";
 
CenterDataVector::usage=
"CenterDataVector[x] the vector is rescaled to [0,1].";
 
SphereDataMatrix::usage=
"SphereDataMatrix[x] each column is transformed by mean correction
and dividing by the standard deviation (divisor 1/n is used for the
sample variance).";
 
SphereDataVector::usage=
"SphereDataVector[x] the vector is transformed by mean correction
and dividing by the standard deviation (divisor 1/n is used for the
sample variance).";

RemoveColumns /: RemoveColumns::usage = 
     "RemoveColumns[x] removes columns specified by i from data matrix x.";

FindRedundantDimensions /: FindRedundantDimensions::usage = "FindRedundantDim\
ensions[x] takes a data matrix with n rows and p columns and finds any \
columns containing a constant value.";
 
Randu /: Randu::usage = "Randu[] generates one uniform (0,1) pseudorandom \
number using the RANDU generator. Randu[n] generates n uniform(0,1) numbers. \
If no seed is supplied, the initial seed is set to: RanduSeed=1345577999.  \
Alternatively you can initial choose any integer value for RanduSeed.";

(*----------------------------------------------------------------------*)
(*Parallel Coordinates Plot*)

Options[ParallelCoordinatesPlot] = {ParallelLines -> True, 
     ShowPoints -> False, ShowLines -> True, ColorLines -> False, 
     CommonScale -> False, AxesLabels -> False};
 
ParallelCoordinatesPlot /: ParallelCoordinatesPlot::redundantDimensions = "Co\
lumn(s) `1` contain constant values and supply no information -- they are \
redundant.";
 
ParallelCoordinatesPlot /: ParallelCoordinatesPlot::removed = 
     "Redundant columns removed before centering and standardizing to [0,1].";
 
ParallelCoordinatesPlot /: ParallelCoordinatesPlot::usage = "ParallelCoordina\
tesPlot[x] produces a parallel coordinates plot for data matrix x. The \
default values of the options are: \
ParallelLines\[Rule]True,ShowPoints\[Rule]False,ShowLines\[Rule]True,ColorLin\
es\[Rule]False,CommonScale\[Rule]False,AxesLabels\[Rule]False.  To omit \
showing the parallel lines, set ParallelLines->False.  To omit joining the \
points with lines, set ShowLines->False.  In order to show the points, set \
ShowPoints->True.  To color the lines joining the points set ColorLines equal \
to a vector of colors of length n, where n is the number of rows of x.  To \
use the same scaling on each axis, set CommonScale->True.  To label the axes, \
set AxesLabels equal to a vector of character strings of length p, where p is \
the number of columns of x.  To select the observations corresponding to the \
k-most extreme values in column j, set SelectExtremes->{k,j}.  When \
SelectExtremes is used the observations corresponding to high values are in \
red and the low ones are in blue.  Also see LineLegend for a creating a \
separate graphic for the coloring used when ColorLines is set.";
 

LineLegend::usage=
"LineLegend[colors, titles] produces horizontal lines specificed by 
colors and labelled with titles. Suitable for creating legends.";

(*----------------------------------------------------------------------*)
(*input-output*)

Options[ReadData] = {Silent -> False};

ReadData/: ReadData::usage = "ReadData[file] file is a character string which \
specifies a file. ReadData[file] inputs a  data file. \
Initial lines in the file which begin with the character # are assumed to contain \
comments and documentation. The remainder of the file contains the data values \
\!\(z\_1\),\!\(z\_2\), ... which are read by rows.  \
The output is a vector containing all the data. \
There is one option which has default value, Silent->False. 
When Silent->False, the title is displayed as well as the number of data values input. \
When Silent->True, the title is not displayed nor the number of data values input.\
If the file contains p columns and n rows, you can obtain the n-by-p matrix by using
Partition.  For example, Partition[ReadData[file],p] returns an n-by-p matrix
assuming that there are n*p data values present in the file.";
 
Options[WriteData] = {Title -> False};

WriteData/: WriteData::usage = "WriteData[filename, x] outputs the numeric list \
x to filename. If x is not a vector it is flattened and then output as one
vector.  WriteData takes an option argument, Title. If Title is set a \
character string then this string is first on  the output file. The \
default setting for Title is False.";


(*==========================================================================*)
Begin["`Private`"]

(* utility *)
NumericVectorQ[x_] := VectorQ[x] && And @@ NumericQ /@ x;

NumberVectorQ[x_] := VectorQ[x] && And @@ NumberQ /@ x;

NumberListQ[x_] := ListQ[x] && And @@ NumericQ /@ Flatten[x];

PowerTransformation[z_, lambda_] := 
    If[Abs[lambda] < 10^(-7), Log[z], z^lambda];

RescaleDataMatrix[x_] := StandardizeDataMatrix[SphereDataMatrix[x]];
 
StandardizeDataMatrix[x_] := Transpose[StandardizeDataVector /@ Transpose[x]];
 
StandardizeDataVector[x_] := (x - Min[x])/(Max[x] - Min[x]);
 
CenterDataMatrix[x_] := Transpose[CenterDataVector /@ Transpose[x]];
 
CenterDataVector[x_] := x - Plus @@ x/Length[x];

SphereDataMatrix[x_] := Transpose[SphereDataVector /@ Transpose[x]];
 
SphereDataVector[x_] := With[{mean = Plus @@ x/Length[x]}, 
     (x - mean)/Sqrt[Plus @@ ((x - mean)^2)/Length[x]]];

FindRedundantDimensions[x_] := BooleanSelect[Range[Last[Dimensions[x]]], 
     (And @@ MapThread[Equal, Transpose[Partition[#1, 2, 1]]] & ) /@ 
      Transpose[x]];
 
RemoveColumns[x_, i_] := Transpose[Transpose[x][[
      Select[Range[Last[Dimensions[x]]],  !MemberQ[i, #1] & ]]]];

Randu[] := Module[{}, RanduSeed = If[NumberQ[RanduSeed], RanduSeed, 
        RanduSeed = 1345577999]; RanduSeed = Mod[(3 + 2^16)*RanduSeed, 2^31]; 
      RanduSeed/2.^31];
 
Randu[n_] := Table[Randu[], {n}];

(*----------------------------------------------------------------------*)
(*ParallelCoordinatesPlot*)

ParallelCoordinatesPlot[x_, (opts___)?OptionQ] := 
    Module[{p, n, g = {}, prllines, colors, spoints, slines, cscale, 
        LineColor, extremes, y = x, i, ind, je, ke, lo, hi, caseslow, caseshi,
         ye, adjx = 0.05}, 
      alabels = AxesLabels /. {opts} /. Options[ParallelCoordinatesPlot];
      colors = ColorLines /. {opts} /. Options[ParallelCoordinatesPlot];
      cscale = CommonScale /. {opts} /. Options[ParallelCoordinatesPlot];
      extremes = SelectExtremes /. {opts} /. Options[ParallelCoordinatesPlot];
       ind = FindRedundantDimensions[x];
      If[Length[ind] > 0, 
        Message[ParallelCoordinatesPlot::redundantDimensions, ind]];
      y = 
        If[! cscale && Length[ind] > 0, 
          Message[ParallelCoordinatesPlot::removed, {}]; 
          RemoveColumns[y, ind], y]; {n, p} = Dimensions[y];
      y = If[cscale, (y - Min[y])/(Max[y] - Min[y]), RescaleDataMatrix[y]];
      LineColor = 
        If[Length[colors] == 0, Table[RGBColor[0, 0, 0], {n}], colors]; 
      prllines = ParallelLines /. {opts} /. Options[ParallelCoordinatesPlot]; 
      spoints = ShowPoints /. {opts} /. Options[ParallelCoordinatesPlot]; 
      slines = ShowLines /. {opts} /. Options[ParallelCoordinatesPlot];
      g = 
        If[prllines, {GrayLevel[0.7], Thickness[0.02], 
            Line /@ Partition[
                Partition[
                  Flatten[{Transpose[{Range[p], Array[0 &, p], Range[p], 
                          Array[1 &, p]}]}], 2], 2]}, g];
      y = If[Length[extremes] == 0, y, {je, ke} = extremes;
          {lo, hi} = Sort[ye = Transpose[y][[je]]][[{ke + 1, n - ke}]];
          caseslow = BooleanSelect[Range[n], (#1 < lo &) /@ ye];
          caseshi = BooleanSelect[Range[n], (#1 > hi &) /@ ye];
          y[[Join[caseslow, caseshi]]]]; 
      LineColor = 
        If[Length[extremes] == 0, LineColor, 
          Flatten[Transpose[
              Table[{RGBColor[0, 0, 1], RGBColor[1, 0, 0]}, {ke}]]]]; {n, p} =
         Dimensions[y];
      g = 
        If[slines, {g, 
            Table[{LineColor[[i]], Line[Transpose[{Range[p], y[[i]]}]]}, {i, 
                n}]}, g];
      g = 
        If[spoints, {g, 
            Table[Point /@ Transpose[{Range[p], y[[i]]}], {i, n}]}, g]; 
      g = If[Length[alabels] > 0, {g, 
            MapThread[
              Text[#1, #2, {0, 0}] &, {alabels, 
                Transpose[{Range[p], Array[-0.1 &, p]}]}], {RGBColor[1, 1, 1],
               Line[{{1 - adjx*p, 0}, {(1. + adjx)*p, 0}}]}}, g];
      Show[Graphics[g], TextStyle -> {FontWeight -> "Bold", FontSize -> 8}]];

LineLegend[colors_, titles_] := 
    Module[{g1, g2, K = Length[colors]}, 
      g1 = {{RGBColor[1, 1, 1], Line[{{0, 0.8}, {3, 0.8}}], 
            Line[{{0, K*1.1}, {3, K*1.1}}]}, 
          Transpose[{colors, 
              Line /@ Transpose[
                  Transpose /@ {{Array[0 &, K], Range[K]}, {Array[0.8 &, K], 
                        Range[K]}}]}]};
      g2 = 
        MapThread[
          Text[#1, #2, {-1, 0}] &, {titles, 
            Transpose[{Array[1 &, K], Range[K]}]}];
      Show[Graphics[{AbsoluteThickness[2], g1, g2}], ImageSize -> 144, 
        TextStyle -> {FontWeight -> "Bold", FontSize -> 12}]];

(*----------------------------------------------------------------------*)
(*input-output*)

ReadData[(datafile_)?StringQ, (opts___)?OptionQ] := 
    Module[{continuetitle = True, z = {}, i, c, title, zt, NoTitle = True}, 
     sQ = Silent /. {opts} /. Options[ReadData]; OpenRead[datafile]; 
      title = ""; While[continuetitle, i = StreamPosition[datafile]; 
        c = ToString[Read[datafile, Character]]; 
        c = Which[LetterQ[c], c, c == "#", "", True, 0]; 
        If[c != 0, NoTitle = False]; If[c == 0, continuetitle = False, 
         title = StringJoin[StringJoin[title, "\n"], Read[datafile, String, 
            RecordSeparators -> "\n"]]]]; SetStreamPosition[datafile, i]; 
      While[ToString[zt = Read[datafile, Real]] != "EndOfFile", 
       AppendTo[z, zt]]; Close[datafile]; 
      If[ !sQ, If[ !NoTitle, Print[StringDrop[title, 1]]]; 
               Print["Number of Values Input = ", Length[z]]
        ]; 
     z];
 
WriteData[(filename_)?StringQ, (x_)?NumberListQ, (opts___)?OptionQ] := 
    Module[{x2 = N[x//Flatten, 10], title}, 
     title = Title /. {opts} /. Options[WriteData]; 
     title = If[StringQ[title], 
                StringReplace[StringJoin["#",title], "\n" -> "\n#"],
                title];
     OpenWrite[filename]; 
     If[StringQ[title], 
        WriteString[filename, StringJoin[title, "\n", 
          StringReplace[StringDrop[StringDrop[ToString[x2], 1], -1], 
          "," -> "\n"]]], 
        WriteString[filename, StringReplace[StringDrop[StringDrop[
          ToString[x2], 1], -1], "," -> "\n"]]]; 
     Close[filename]];


  End[];

EndPackage[]




