(*:Version: Mathematica 4.0 *)

(*:Name: Statistics`Common`MultivariateCommon` *)

(*:Context: Statistics`Common`MultivariateCommon` *)

(*:Title: Symbols for Multivariate Statistics *)

(*:Author:
  E.C. Martin (Wolfram Research), 1994
*)

(*:Copyright: Copyright 1994-1999, Wolfram Research, Inc. *)

(*:History:
   Original package, 1994, ECM.
   Graphics rules made more efficient, 1998, ECM.
*)

(*:Reference: Usage messages only. *)

(*:Summary:
This package defines symbols in a common context for use in multivariate
statistical packages.
*)

(*:Requirements: No special system requirements. *)

(*:Warning: None. *)

(*:Sources: Basic statistics texts. *)

BeginPackage["Statistics`Common`MultivariateCommon`"]
 
CovarianceMatrix

CorrelationMatrix

(* geometric primitives *)

Ellipsoid::usage =
"Ellipsoid[{x1, ..., xp}, {r1, ..., rp}, {d1, ..., dp}] represents a
p-dimensional ellipsoid centered at the point {x1, ..., xp}, where the ith
semi-axis has radius ri and lies in direction di.  Ellipsoid[{x1, ..., xp},
{r1, ..., rp}, IdentityMatrix[p]] simplifies to Ellipsoid[{x1, ..., xp},
{r1, ..., rp}], an ellipsoid aligned with the coordinate axes."

Polytope::usage =
"Polytope[{{x11, ..., x1p}, ..., {xn1, ..., xnp}}, connectivity] represents
a p-dimensional polytope, with n vertices {x11, ..., x1p}, ..., {xn1, ..., xnp},
where the connections between vertices is specified by connectivity."

Begin["`Private`"]

(* Eliminate third argument specifying
         direction if the Ellipsoid is aligned with coordinate axes. *)
Ellipsoid[mu_?VectorQ, r_?VectorQ, dir_?MatrixQ] :=
  Module[{rdir = Transpose[{r, dir}], p = Length[mu],
          sortedR, sortedDIR, newEllipsoid},
           newEllipsoid /;   (
                      {sortedR, sortedDIR} =
                        Transpose[Sort[rdir, Order[#1[[2]], #2[[2]]]==-1&]];
                      If[sortedDIR[[1, 1]] != 0,
                         sortedDIR /= sortedDIR[[1, 1]] ];
                      If[Apply[And, Map[#==0&, Flatten[
                                sortedDIR - IdentityMatrix[p] ]]],
                         newEllipsoid = Ellipsoid[mu, sortedR]; True,
                         False
                      ])
  ] /; Apply[Equal, Join[Dimensions[dir], {Length[mu], Length[r]}]]

(* For p=2, Graphics[Ellipsoid[]] evaluates.
   For p=3, Graphics3D[Ellipsoid[]] evaluates.
   For p=1 or p>3, no graphics is defined.
*)

(* p=2 *)
Ellipsoid/: Graphics[Ellipsoid[mu_, r_, dir_?MatrixQ], opts___Rule] :=
  graphicsEllipsoid[mu, r, dir[[1]], opts] /; (Length[mu] == Length[r] ==
                                        Dimensions[dir][[2]] == 2) &&
                                        1 <= Dimensions[dir][[1]] <= 2

(* For p=2, permit the direction to be specified as {a, b}, in addition to
        {{a, b}} and {{a, b}, {c, d}}. *)
Ellipsoid/: Graphics[Ellipsoid[mu_, r_, d1_?VectorQ], opts___Rule] :=
  graphicsEllipsoid[mu, r, d1, opts] /; Length[mu] == Length[r] ==
                                        Length[d1] == 2


(* For p=2, use Circle instead of ParametricPlot when ellipse is not tilted. *)
Ellipsoid/: Graphics[Ellipsoid[mu_, r_], opts___Rule] :=
  Graphics[Circle[mu, r], opts] /; Length[mu]==Length[r]==2


(* p=3 *)
Ellipsoid/: Graphics3D[Ellipsoid[mu_, r_, dir_?MatrixQ], opts___Rule] :=
  Module[{plotpoints},
     plotpoints = PlotPoints /. {opts} /. {PlotPoints -> Automatic};
     (* assume all direction vectors have unit length *)
     graphicsPrimitives = First[ParametricPlot3D[mu +
        Transpose[dir] . (r {Cos[t] Cos[u], Sin[t] Cos[u], Sin[u]}),
         {t, 0, 2Pi}, {u, -Pi/2, Pi/2},
         DisplayFunction->Identity, PlotPoints->plotpoints, Compiled->False]];
     Graphics3D[graphicsPrimitives, opts]
  ] /; (Length[mu] == Length[r] ==
                        Dimensions[dir][[1]] == Dimensions[dir][[2]] == 3)

(* Graphics[Ellipsoid[]] code by Jeff Adams, 10/93 *)
graphicsEllipsoid[mu_, r_, d1_, opts___Rule] :=
  Module[{cos, sin, plotpoints, graphicsPrimitives},
     plotpoints = PlotPoints /. {opts} /. {PlotPoints -> 30};
     (* assume d1 has unit length *)
     {cos, sin} = d1;
     graphicsPrimitives = First[ParametricPlot[mu +
         {{cos, -sin}, {sin, cos}} .
         (r {Cos[theta], Sin[theta]}), {theta, -Pi, Pi},
         DisplayFunction->Identity, PlotPoints->plotpoints, Compiled->False]];
     Graphics[graphicsPrimitives, opts]
  ]

(* ===================================================================== *)

Format[Polytope[vertices_, connectivity_List]] :=
        Polytope[vertices, "-Connectivity-"]

(* For p=2, Graphics[Polytope[]] evaluates.
   For p=3, Graphics3D[Polytope[]] does not evaluate because 3D convex hulls
        are not yet implemented.
   For p=1 or p>3, no graphics is defined.
*)

(* p=2 *)
Polytope/: Graphics[Polytope[vertices_?MatrixQ, connectivity_], opts___] :=
        Graphics[Line[Append[vertices, First[vertices]]], opts] /;
                        Dimensions[vertices][[2]]==2


(* ================ rules permitting graphics directives =================== *)

Unprotect[Graphics]
Graphics[{d1___, Ellipsoid[e___], d2___}, opts___Rule] :=
        With[{primitive1 = Graphics[Ellipsoid[e]][[1]]},
                Graphics[Join[{d1}, {primitive1}, {d2}], opts]
        ]

Graphics[{d1___, Polytope[p___], d2___}, opts___Rule] :=
        With[{primitive1 = Graphics[Polytope[p]][[1]]},
                Graphics[Join[{d1}, {primitive1}, {d2}], opts]
        ]

Graphics[{g1___, {directives___, Ellipsoid[e___]}, g2___}, opts___Rule] :=
        With[{primitive1 = Graphics[Ellipsoid[e]][[1]]},
                Graphics[{g1, Join[{directives}, {primitive1}], g2}]
        ]

Graphics[{g1___, {directives___, Polytope[p___]}, g2___}, opts___Rule] :=
        With[{primitive1 = Graphics[Polytope[p]][[1]]},
                Graphics[{g1, Join[{directives}, {primitive1}], g2}]
        ]
Protect[Graphics]

Unprotect[Graphics3D]
Graphics3D[{d1___, Ellipsoid[e___], d2___}, opts___Rule] :=
        With[{primitive1 = Graphics3D[Ellipsoid[e]][[1]]},
                Graphics3D[Join[{d1}, {primitive1}, {d2}], opts]
        ]

Graphics3D[{g1___, {directives___, Ellipsoid[e___]}, g2___}, opts___Rule] :=
        With[{primitive1 = Graphics3D[Ellipsoid[e]][[1]]},
                Graphics3D[{g1, Join[{directives}, {primitive1}], g2}]
        ]
Protect[Graphics3D]


End[]

EndPackage[]
