info     0
) Graphical macros useful in Multivariate Analysis
)
) Plots of multivariate data (see also plotmatrix() in Graphics.mac)
)  andrewsplot()      Make Andrews plot
)  parcoordplot()     Make parallel coordinate plot
)  faces()            Make Chernoff-like faces
)
) Diagnostic plots
)  chiqqplot()        Plots of quadratic form in data or residuals
)                     against chi-squared probability points
)
) Cluster analysis-related plots
)  cluscritplot()     Plot of criterion or changes to help choose the
)                     number of clusters
)
) Supporting macros
)  interpspline()     Macro to compute interpolating spline; used by
)                     faces;  it or an equivalent will probably be
)                     moved to Math.mac someday
) Help is at the end of the file
))040910 Version
))040918 Modified so that faces(...,oldversion:T) will work on a version
))       of MacAnova that doesn't implement str$comp <- rhs
===> andrewsplot <===
andrewsplot macro dollars
) Macro to make parallel coordinate plots of multivariate data
) Usages:
)  andrewsplot(x [,rescale:T or stand:T] [,npoints:np] [,graph keys])
)  andrewsplot(x, groups [,linetypes:ltype] [,rescale:T or stand:T]\
)           [,npoints:np] [,graph keywords])
)   x          REAL N by p data matrix containing no MISSING values
)   groups     Vector of positive integers defining groups or samples,
)              with length(groups) = N
)   ltype      Length g vector of positive integers specifying line
)              types, where g = number of distinct values in groups;
)              default is ltype = run(g)
)   np         positive integer > 1, number of points at which curves
)              are plotted; default np = 35
)   rescale:T  Rescale as x -> (x - min(x))/(max(x)-min(x)
)   stand:T    Standardize: x -> (x - mean(x))/stddev(x)
)   graph keys add:T, show:F, window:m and any graphics keyword phrases
)              that can be used with showplot()
)
)  If y is x or a rescaled or standardized version of x, functions
)  f_i(t) are computed as trigonimetric series with coefficients from
)  row y[i,] and plotted as a line graphs for i = 1,...,nrows(x).
)
)  Specifically the function for case i is
)
)   f_i(t) = y[i,1]/sqrt(2) + y[i,2]*cos(t) + y[i,3]*sin(t) +
)       y[i,4]*cos(2*t) ... ,  -.5 <= t <= +.5)  
)
)  The last term is y[i,p]*sin(m*t) or y[i,p]*cos(m*t) depending on
)  whether p = ncols(x) = 2*m is even or p = 2*m+1 is odd
)
)  The arguments to cos() and sin() are in cycles (1 = 2*PI radians = 360 degrees).
)
)  When groups is an argument, different line types are used for
)  different groups.
)
))031206 Written by C. Bingham kb@umn.edu
))031207 Modified so it doesn't go in a filled window unless add:T or
))       window:m are arguments
#$S(x [,rescale:F or stand:T] [,graph keywords])
#$S(x, groups [,linetypes:ltype] [,rescale:F or stand:T] [,graph keywords])
@y <- matrix(argvalue($1,"x","nonmissing real matrix"))
@p <- ncols(@y)
@n <- nrows(@y)
@groups <- if ($v > 1){
	argvalue($2,"groups","positive integer vector")
} else {
	rep(1,@n)
}
if (@n != nrows(@groups)){
	error("length($2) != nrows($1)")
}
@g <- max(@groups)

@rescale <- keyvalue($K,"resc*","TF",default:F)
@stand <- keyvalue($K,"stand*","TF",default:F)
if (@rescale && @stand) {
	error("illegal to use 'rescale:T' with 'stand:T'")
}
@npoints <- keyvalue($K,"npoint*","positive count",default:35)
if (@npoints < 2) {
	error("value for 'npoints' must be at least 2")
}
@add <- keyvalue($K,"add","TF",default:F)
@show <- keyvalue($K,"show","TF")
@window <- keyvalue($K,"win*","count")
if (isnull(@show)) {
	@show <- T
} elseif (!@show && !isnull(@window)) {
	error("'window' used with 'show:F'")
}

@linetypes <- keyvalue($K,"linetype*","positive integer vector",\
	default:run(@g))

if (length(@linetypes) < @g){
	error("Fewer line types than groups")
}

if (@rescale) {
	@y <- hconcat((@y - min(@y))/(max(@y)-min(@y)))
} elseif (@stand) {
	if (!ismacro(standardize)) {
		getmacros(standardize,silent:T)
	}
	@y <- standardize(@y)
}
@u <- run(-.525,.525,1.05/(@npoints-1))
@U <- rep(1/sqrt(2),@npoints) + rep(0,@p)'

 # @U is hconcat(rep(1/sqrt(2),npoints),cos(t),sin(2t),cos(2*t)...
if (@p > 1) {
	for (@i,2,@p,2) {
		@j <- @i/2
		@U[,@j+1] <- sin(@j*@u,cycles:T)
		if (@j+1 < @p) {
			@U[,@j+2] <- cos(@j*@u,cycles:T)
		}
	}
	delete(@i,@j)
}

if (@n > 1) {
	# omit last case when plotting to a specific window or not displaying 
	@n1 <- @n - (!isnull(@window) || !@show)
	for (@i, 1, @n1) {
		lineplot(@u,@U %C% @y[@i,],linetype:@linetypes[@groups[@i]],show:F,\
			     add:@add || @i > 1)
	}
	delete(@i,@n1)
}

@title <- if (@stand) {
	"Standardized Andrews plot"
} elseif (@rescale) {
	"Rescaled Andrews plot"
} else {
	"Andrews plot"
}
if (@g > 1) {
	@title <- paste(@title,"split by $2")
}

if (!isnull(@window) || !@show) { # plot last case
	lineplot(@u,@U %C% @y[@n,],linetype:@linetypes[@groups[@n]],\
		   $K, xticks:run(-.5,.5,.1),xmin:-.5,xmax:.5, ymin:?,ymax:?,\
		   xaxis:F,yaxis:F,title:@title,xlab:"t",add:@add || @n > 1)
} else { #display graph
	showplot(show:T,$K, xticks:run(-.5,.5,.1),xmin:-.5,xmax:.5, ymin:?,ymax:?,\
		   xaxis:F,yaxis:F,title:@title,xlab:"t")
}
delete(@u,@U,@y,@groups,@linetypes,@stand,@rescale,@n,@g,@p,@title,\
	   @add,@show,@window)
%andrewsplot%

===> parcoordplot <===
parcoordplot macro dollars
) Macro to make parallel coordinate plots of multivariate data
) Usages:
)  parcoordplot(x [,rescale:F] [,graph keys])
)  parcoordplot(x, stand:T [,graph keys])
)  parcoordplot(x, groups [,linetypes:ltype] [,rescale:F] [,graph keys])
)  parcoordplot(x, groups, stand:T [,linetypes:ltype] [,rescale:F]\
)           [,graph keys])
)   x          REAL data matrix with at least 2 columns contining
)              non-MISSING values
)   groups     Vector of positive integers defining groups or samples,
)              with length(groups) = nrows(x)
)   ltype      Length g vector of positive integers specifying line
)              types, where g = number of distinct values in groups;
)              default is ltype = run(g)
)   rescale:F  Don't rescale as x -> (x - min(x))/(max(x)-min(x);
)              default is to rescale
)   stand:T    Standardize: x -> (x - mean(x))/stddev(x) instead of
)              rescaling
)   graph keys add:T, show:F, window:m or any graphics keyword phrases
)              that can be used with showplot()
)
)  If y is the rescaled, standardized or unchanged version of x,
)  each row y[i,] is plotted as a line graph vs 1, 2, ..., ncols(x)
)
)  When groups is an argument, different line types are used for
)  different groups.
)
))031206 Written by C. Bingham kb@umn.edu
))031207 Modified so it doesn't go in a filled window unless add:T or
))       window:m are arguments
#$S(x [,rescale:F or stand:T] [,graph keywords])
#$S(x, groups [,linetypes:ltype] [,rescale:F or stand:T] [,graph keywords])
@y <- matrix(argvalue($1,"x","real matrix"))
@p <- ncols(@y)
@n <- nrows(@y)
@groups <- if ($v > 1){
	argvalue($2,"groups","positive integer vector")
} else {
	rep(1,@n)
}
if (@n != nrows(@groups)){
	error("length($2) != nrows($1)")
}
@rescale <- keyvalue($K,"rescale","TF")
@stand <- keyvalue($K,"stand*","TF")

if (isnull(@rescale) && isnull(@stand)) {
	@rescale <- T
	@stand <- F
} elseif (isnull(@stand)) {
	@stand <- F
} elseif (isnull(@rescale)) {
	@rescale <- !@stand
} elseif (@rescale && @stand) {
	error("illegal to use 'rescale:T' with 'stand:T'")
}
@add <- keyvalue($K,"add","TF",default:F)
@show <- keyvalue($K,"show","TF")
@window <- keyvalue($K,"win*","count")
if (isnull(@show)) {
	@show <- T
} elseif (!@show && !isnull(@window)) {
	error("'window' used with 'show:F'")
}

@g <- 	max(@groups)
@gg <- length(unique(@groups))

@linetypes <- keyvalue($K,"linetype*","positive integer vector",\
	default:run(@g))

if (length(@linetypes) < @g){
	error("Fewer line types than groups")
}
if (sum(vector(sum(!ismissing(@y))) > 0) < 2) {
	error("You need at least 2 columns with non-MISSING data")
}

if (haslabels(@y)){
	@xticklabs <- getlabels(@y,2)
	@xlab <- "Variable"
} else {
	@xticklabs <- paste(run(@p),multiline:T)
	@xlab <- "Variable Number"
}
@J <- run(1,@p,ceiling(@p/8))
@xticks <- run(@p)[@J]
@xticklabs <- @xticklabs[@J]

@oldwarn <- getoptions(warnings:T)
setoptions(warnings:F) # suppress warnings about MISSING values
@ymin <- min(@y)
@ymax <- max(@y)
if (@rescale) {
	@y <- hconcat((@y - @ymin)/(@ymax-@ymin))
} elseif (@stand) {
	if (!ismacro(standardize)) {
		getmacros(standardize,silent:T)
	}
	@y <- standardize(@y)
}
@ymin <- min(vector(@y))
@ymax <- max(vector(@y))
setoptions(warnings:@oldwarn) # restore option value

@x <- rep(run(@p+1), @n)
@y <- vector(hconcat(@y,rep(?,@n))')
if (@gg > 1 || @show || !isnull(@wind)){
	# omit last group when plotting to a specific window or not displaying 
	@g1 <- @g - (!isnull(@window) || !@show)
	for(@i,1,@g1){
		@J <- @groups == @i
		if (sum(@J) > 0){
			@J <- 1 + vector(((@p+1)*run(0,@n-1)[@J]+run(0,@p)')')
			lineplot(@x[@J],@y[@J],linetype:@linetypes[@i],show:F,\
				   add:@add || @i>1)
		}
	}
	delete(@i,@g1)
}

@title <- if (@stand) {
	"Standardized parallel coordinates plot"
} elseif (@rescale) {
	"Rescaled parallel coordinates plot"
} else {
	"Parallel coordinates plot"
}
if (@gg > 1) {
	@title <- paste(@title,"split by $2")
}

if (!@show || !isnull(@window)) {
	# plot last group
	@J <- 1 + vector(((@p+1)*run(0,@n-1)[@groups == @g]+run(0,@p)')')
	lineplot(@x[@J],@y[@J],linetype:@linetypes[@g],add:@add || @gg > 1,$K,\
		   xticks:@xticks,xticklabs:@xticklabs,xmin:1,xmax:@p,\
		   ymin:1.05*@ymin-.05*@ymax,ymax:1.05*@ymax-.05*@ymin,\
		   xaxis:!(@stand||@rescale),title:@title,xlab:@xlab)
} else {
	showplot(show:T,$K,xticks:@xticks,xticklabs:@xticklabs,xmin:1,xmax:@p,\
		   ymin:1.05*@ymin-.05*@ymax,ymax:1.05*@ymax-.05*@ymin,\
		   xaxis:!(@stand||@rescale),title:@title,xlab:@xlab)
}

delete(@x,@J,@ymin,@ymax,@title,@xlab,@xticks,@xticklabs,@add)
delete(@y,@groups,@n,@g,@gg,@p,@stand,@rescale,@linetypes)
%parcoordplot%


===> interpspline <===
interpspline    macro   dollars
) macro to compute spline interpolant
) Usage
)  result <- interpspline(a, y [, npoints:m] [,plot:T])
)   a         nonMISSING REAL vector
)   y         nonMISSING REAL vector
)   m         integer > 0, default:50
)   plot:T    make line plot
)   result    hconcat(x, Sx), x = a[1] + (a[2]-a[1])*run(0,m-1)/(m-1),
)             Sx = REAL vector of length m containing spline values
)
) The algorithm is from Ralson, A and Rabinowitz, P (1985), A First
) Course in Numerical Analysis, McGraw-Hill, pp 76ff
)
)) 040112 Translated from R code of Dr. Hans Peter Wolf
)) (pwolf@wiwi.uni-bielefeld.de), downloaded 10 Jan 2004 from
)) http://www.wiwi.uni-bielefeld.de/~wolf/software/R-wtools/faces/
#$S(a, y [, npoints:m] [,plot:T]), a, y REAL vectors, m>0 integer
@a <- argvalue($1,"a","nonmissing real vector")
@y <- argvalue($2,"y","nonmissing real vector")
@m <- keyvalue($K,"npoint*","positive count",default:50)
@plot <- keyvalue($K,"plot","TF",default:F)
@n <- length(@a)
if (length(@y) != @n){
	error("Arguments 1 and 2 have different lengths")
}
@h <- movavg(1,@a)[-1]
@dy <- movavg(1,@y)[-1]
@sigma <- @dy/@h
@hh <- @h[-1] + @h[-(@n-1)]  # @h[1]+2*(@h[2] + ... + @h[@n-3]) + @h[@n-2]
@lambda <- @h[-1]/@hh
@mu <- 1 - @lambda
@d <- 6*movavg(1,@sigma)[-1]/@hh

 # construct n-2 by n-2 tri-diagonal matrix
@tri_mat <- 2*dmat(@n-2,1)
@tri_mat[run(1,@n-3) + vector(1,0)'] <- @mu[-1]          # below diagonal
@tri_mat[run(1,@n-3) + vector(0,1)'] <- @lambda[-(@n-2)] # above diagonal

@M <- vector(0,solve(@tri_mat, @d),0)

@x <- run(@a[1],@a[@n],(@a[@n]-@a[1])/(@m-1))
@ans_kl <- bin(@x - 1e-12*(@x - @a[@n]),@a)$counts

@i <- rep(run(@n-1),@ans_kl) + 1
@i1 <- @i - 1
@S_x <- @M[@i1]*(@a[@i] - @x)^3 / (6*@h[@i1])  + \
		@M[@i]*(@x - @a[@i1])^3 / (6*@h[@i1])  + \
		(@y[@i1] - @M[@i1]*@h[@i1]^2/6) * (@a[@i] - @x)/@h[@i1] + \
		(@y[@i] - @M[@i]*@h[@i1]^2/6) * (@x - @a[@i1]) /@h[@i1]

if (@plot){
	lineplot(@x,@S_x, show:F)
	addpoints(@a,@y,$K)
}
delete(@plot,@a,@y,@i1,@i,@M,@h,@hh,@tri_mat,@d,@mu,@lambda)
hconcat(delete(@x,return:T),delete(@S_x,return:T),\
		labels:structure("@",vector("x","spline(x)")))
%interpspline%


===> faces <===
faces macro dollars
) Usage:
)  faces(Y [, whichvars:varnos] [,byrow:T] [,mods:Mods]\
)          [,scale:F] [,min:Max] [,max:Max]\
)          [,keep:T] [,draw:T] [,nrows:nr] [,ncols:nc] [,panel:F] \ 
)          [,labels:facelabs] [,labelplace:lp]\
)          [,gridlines:T] [,npoints:npts] [graph keywords])
)   Y          n by m REAL data matrix with no missing values;
)              rows are cases, columns are variables
)   varnos     vector of positive integers <= m; faces will be created
)              from rows of Y[,varnos], default is run(m)
)   byrow:T    means columns are cases and rows are variables
)   Mods       structure with positive scalar components with
)              names chosen from 'facesize' [.2], 'faceshape' [.2],
)              'mouthsize' [.7], 'smile' [20], 'eyesize' [.7],
)              'hairsize' [.2], 'hairstyle' [50], 'noseheight' [.7],
)              'nosewidth' [1], or 'earsize' [.7] used to replace
)              default values (in [...]) used to translate variables
)              into feature characteristics
)   Min, Max   REAL column or row vectors of length m or scalars, to be
)              used in scaling in places of min(Y) and/or max(Y)
)   scale:F    Don't rescale variables to [-1,1] interval; values < -1
)              are rounded to -1 and values > +1 are rounded to +1;
)              default is to scale to (Y - Min)/(Max - Min)
)   keep:T     Return coordinates of faces in structure of matrices with
)              two columns, one matrix for each face, instead of drawing
)              faces
)   draw:T     Draw faces, even with keep:T
)   panel:F    Don't make panel graph; make separate plot for each face
)   nr, nc     positive integers, numbers of rows and columns in
)              array of faces; ignored with 'panel:F'
)   gridline:T Draw framing faces; ignored with 'panel:F'
)   facelabs   length n CHARACTER vector of labels for faces; default
)              is "Case xx" with panel:F and "#" (use case numbers)
)              without panel:F.  With keep:T, the component names of the
)              result will be facelabs (default "Case_xx") which should
)              contain no spaces or '$'
)   lp         CHARACTER scalar specifying position of labels; default
)              is "midleft" or "ml"; ignored with panel:F.  The first
)              letter defines vertical position, "b" (bottom), "m" (mid),
)              or "t".  A later "l" (left), "c" (center) or "r" (right)
)              defines a horizontal position.  "mc" is not a good idea.
)   npts       Integer > 0, number of points in interpolating spline used
)              for each curve drawn
)
) You can use any graphics keywords that can be used with panelplot()
)
) This requires a version of panelplot() dated 040112 or later
)           Features parameters of this implementation           Default
)  1-height of face   y = y*(1 + facesize*X[1]), all features      [.2]
)  2-width of face    x = x*(1 + facesize*X[2]), all features      [.2]
)  3-shape of face    x = x + (y + 50)*faceshape*sin(-1.5*X[3])    [.2]
)  4-height of mouth  y = m + (y-m)*(1 + mouthsize*X[4]), m=mean   [.7]
)  5-width of mouth   x = x*(1 + mouthsize*X[5])                   [.7]
)  6-curve of smile   x = x + smile*X[6])                          [20]
)  7-height of eyes   y = m + (y-m)*(1 + eyesize*X[7]), m=mean     [.7]
)  8-width of eyes    x = m + (x-m)*(1 + eyesize*X[8]), m=mean     [.7]
)  9-height of hair   y = m + (y-m)*(1 + hairsize*X[9]), m=min     [.2]
) 10-width of hair    x = m + (x-m)*(1 + hairsize*X[10]), m=0      [.2]
) 11-styling of hair  x = x + hairstyle*X[11])                     [50]
) 12-height of nose   y = m + (y-m)*(1 + noseheight*X[12]), m=mean [.7]
) 13-width of nose    x = x*(1 + nosewidth*X[13])                  [1]
) 14-width of ears    x = m + (x-m)*(1 + earsize*X[14]), m=mean    [.7]
) 15-height of ears   y = m + (y-m)*(1 + earsize*X[15]), m=min     [.7]
)
)) 040112 Translated and modified from R code of Dr. Hans Peter Wolf
)) (pwolf@wiwi.uni-bielefeld.de), downloaded 10 Jan 2004 from
)) http://www.wiwi.uni-bielefeld.de/~wolf/software/R-wtools/faces/
)) 040113 added keywords 'mods', 'keep', 'draw', 'labelplace'; corrected
))        some bugs
)) 040114 changed 'reorder' to 'whichvars', requiring only that the
))        values be between 1 and ncols(xy)
)) 040123 Added keywords 'min' and 'max' to modify scaling
)) 040918 added keyword phrase oldversion:T and changed code to allow
))        execution with version that doesn't recognize str$comp <- rhs
#$S(Y [, whichvars:varnos] [,byrow:T] [,mods:Mods]\
#     [,scale:F] [,min:Max] [,max:Max]\
#     [,keep:T] [,draw:T] [,nrows:nr] [,ncols:nc] [,panel:F] \ 
#     [,labels:facelabs] [,labelplace:lp]\
#     [,gridlines:T] [,npoints:npts] [graph keywords])
@xy <- $01
if (isnull(@xy)) { #use artificial example
	@xy <- hconcat(run(3),run(5,3),run(3,5),run(5,7))'
} else {
	@xy <- matrix(argvalue(@xy,"xy","nonmissing real matrix"))
}

@keys <- if ($k > 0) {
	structure($K,notakey:NULL)
} else {
	structure(notakey:NULL)
}

if (keyvalue(@keys,"byrow*","TF",default:F)) {
	 @xy <- @xy'
}

@mm <- ncols(@xy)
@n <- nrows(@xy)

# default modification factors
@mods <- structure(facesize:0.2, faceshape:0.2, mouthsize:0.7,\
				   smile:20, eyesize:0.7, hairsize:0.2, hairstyle:50,\
				   noseheight:0.7, nosewidth:1.0, earsize:0.7)

@whichvars <- keyvalue(@keys,"which*","positive integer vector")
@scale <- keyvalue(@keys,"scale","TF",default:T)
@fill <- keyvalue(@keys,"fill","TF",default:F)
@labels <- keyvalue(@keys,"labels*","character vector")
@nrow <- keyvalue(@keys,"nrow*","positive count")
@ncol <- keyvalue(@keys,"ncol*","positive count")
@panels <- keyvalue(@keys,"panel*","TF",default:T)
@gridlines <- keyvalue(@keys,"grid*","TF",default:F)
@npts <- keyvalue(@keys,"npoint*","positive count")
@labelpl <- keyvalue(@keys,"labelp*","string",default:"midleft")
@keep <- keyvalue(@keys,"keep","TF",default:F)
@draw <- keyvalue(@keys,"draw","TF",default:!@keep)
@min <- keyvalue(@keys,"min","nonmissing real matrix")
@max <- keyvalue(@keys,"max","nonmissing real matrix")
@oldvers <- keyvalue(@keys,"oldv*","TF", default:T)

@min1 <- min(@xy) # row vectors
@max1 <- max(@xy)
if (!isnull(@min)) {
	if (!isvector(@min) && nrows(@min) != 1) {
		error("value for 'min' not row or column vector")
	}
	@min <- matrix(@min,1)

	if (isscalar(@min)) {
		@min <- rep(@min,@mm)'
	} elseif (length(@min) != @mm) {
		error("length of value for 'min' != number of variables")
	}
	if (!@scale) {
		print("WARNING: value for 'min' ignored with 'scale:F'", macroname:T)
	} elseif (min(vector(@min1-@min)) < 0) {
		print("WARNING: at least 1 'min' value > data mininum", macroname:T)
	}
} else {
	@min <- @min1
}

if (!isnull(@max)) {
	if (!isvector(@max) && nrows(@max) != 1) {
		error("value for 'max' not row or column vector")
	}
	@max <- matrix(@max,1)

	if (isscalar(@max)) {
		@max <- rep(@max,@mm)'
	} elseif (length(@max) != @mm) {
		error("length of value for 'max' != number of variables")
	}
	if (!@scale) {
		print("WARNING: value for 'max' ignored with 'scale:F'", macroname:T)
	} elseif (max(vector(@max1-@max)) > 0) {
		print("WARNING: at least 1 'max' value < data maximum", macroname:T)
	}
	if (min(vector(@max-@min)) <= 0) {
		error("at least one value for 'min' >= value for 'max'")
	}
} else {
	@max <- @max1
}

if (isnull(@npts)){
	# use fewer points in spline curves for panel plots of many
	# faces when coordinates are not being kept
	@npts <- if (@keep || !@panels){
		40
	} else {
		min(40,max(10,ceiling(40*(10/@n)))) # 10 <= npts <= 40
	}
}

if (alltrue(@n > 1, (@k <- sum((@max1 == @min1)')) > 0)) {
	print(paste("WARNING:",@k,"variables are constant"), macroname:T)
	delete(@k)
}
delete(@min1,@max1)

if (!@draw){
	@panels <- F
}
if (@panels) {
	# determine vertical and horizontal offsets for labels
	@vert <- if (match("b*",@labelpl,0,exact:F) != 0){
		.07
	} elseif (match("m*",@labelpl,0,exact:F) != 0) {
		.5
	} elseif (match("t*",@labelpl,0,exact:F) != 0) {
		.95
	} else {
		error(paste("\"",@labelpl,"\" doesn't start with t, m or b",sep:""))
	}
	@horiz <- if (match("*l*",@labelpl,0,exact:F) != 0){
		@just <- "l"
		.02
	} elseif (match("*c*",@labelpl,0,exact:F) != 0) {
		@just <- "c"
		.5
	} elseif (match("*r*",@labelpl,0,exact:F) != 0) {
		@just <- "r"
		1
	} else {
		@just <- "c"
		.5
	}

	if (@vert == .5 && @horiz == .5) {
		print("WARNING: you won't like labels in middle of faces")
	}
}

@mods1 <- keyvalue(@keys,"mod*","structure")
if (!isnull(@mods1)) {
	@modnames <- compnames(@mods)
	@modnames1 <- compnames(@mods1)
	for (@i,1,length(@modnames1)){
		@j <- match(paste(@modnames1[@i],"*",sep:""),@modnames,0,exact:F)
		if (@j == 0){
			error(paste("\"",@modnames1[@i],"\" not recognized mod name",\
				  sep:""))
		}
		@mods[@j] <- argvalue(@mods1[@i],paste("value for",@modnames[@j]),\
			"positive number")
	}
	delete(@modnames, @modnames1)
}

delete(@keys, @mods1)

	#4:
@n_char <- 15 # number of characteristics

if (!isnull(@whichvars)) {
	if (max(@whichvars) > @mm) {
		error("maximum 'whichvars' element > number of variables")
	}
	@xy <- @xy[,@whichvars]
	@mm <- ncols(@xy)
	@min <- @min[,@whichvars]
	@max <- @max[,@whichvars]
}
delete(@whichvars)

 #2:
@nc <- if(isnull(@ncol)) {
	ceiling(sqrt(@n))
} else {
	@ncol
}
if(isnull(@nrow)) {
	@nr <- ceiling(@n/@nc)
	if ((@nr - 1)*@nc >= @n) {
		@nr <-- 1
	}
} else {
	@nr <- @nrow
}

if (@nr * @nc < @n) {
	error(paste(@nr,"*",@nc," < ",@n," = number of cases",sep:""))
}

delete(@nrow,@ncol)
 #:2

@defname <- if (@panels) {
	"#"
} else {
	"Case "
}

@xnames <- if (!isnull(@labels)) {
	if (!isscalar(@labels) && length(@labels) != @n) {
		error("length(face_labels) != n")
	}
	@labels
} elseif (haslabels(@xy)) {
	getlabels(@xy,1)
} else {
	@defname
}
if (isscalar(@xnames) && @xnames[1] != "") {
	@xnames <- getlabels(vector(rep(0,@n),labels:@xnames))
}
if (@xnames[1] == "@") {
	@xnames <- getlabels(vector(rep(0,@n),labels:@defname))
}
if (@keep) {
	@defname <- "Case_"
	@xnames1 <- if (!isnull(@labels)) {
		if (!isscalar(@labels) && length(@labels) != @n) {
			error("length(face_labels) != n")
		}
		@labels
	} elseif (haslabels(@xy)) {
		getlabels(@xy,1)
	} else {
		@defname
	}
	if (match("* *",@xnames1, 0, exact:F) != 0 || \
	    match("*$*",@xnames1, 0, exact:F) != 0) {
		print("WARNING: labels can't be used for component names")
		@xnames1 <- @defname
	}
	if (isscalar(@xnames1) && @xnames1[1] != "") {
		@xnames1 <- getlabels(vector(rep(0,@n),labels:@xnames1))
	}
	if (@xnames1[1] == "@") {
		@xnames1 <- getlabels(vector(rep(0,@n),labels:@defname))
	}
}
delete(@labels,@defname)

if(delete(@scale,return:T)){
	@xy <-- @min
	@max <-- @min
	@J <- vector(@max > 0)
	if (sum(@J) > 0) {
		@xy[,@J] <- 2*@xy[,@J]/@max[,@J] - 1
	}
	delete(@J)
} else { #force to [-1,1]
	@xy[vector(@xy) < -1] <- -1
	@xy[vector(@xy) > 1] <- 1
}

if (@mm < @n_char) {
	@xy <- if (!@fill) {
		# extend rows circularly to length n_char
		@J <- rep(run(@mm),ceiling(@n_char/@mm))[run(@n_char)]
		@xy[,delete(@J,return:T)]
	} else {
		#add n_char-mm 0 columns
		hconcat(@xy,matrix(rep(0,@n*(@n_char-@mm)),@n))
	}
} elseif (@mm > @n_char) {
	print(paste("WARNING: last",@mm - @n_char,"variables ignored"))
	# drop mm - n_char columns
	@xy <- @xy[,run(@n_char)]
}


 #:4

#5:

 # eye is 6 point closed curve
 # iris is 4 point closed curve
 # lipso is right half of 9 point closed curve
 # lipsi is right half of 3 point curve which is extended to make
 #   a 5 point curve attached to extreme points of symmetrized lipso
 # nose is right half of 7 point closed curve
 # shape is right half of 17 point closed curve
 # ear is 2 point curve which is extended to make 4 point curve
 #   attached to shape
 # hair is right half of 7 point curve which is extended to make a
 #   9 point curve attached to shape
@face_orig <- structure(\
	eye:hconcat(vector(12,0),vector(19,8),vector(30,8),\
				vector(37,0),vector(30,-8),vector(19,-8),vector(12,0))',\
	iris:hconcat(vector(20,0),vector(24,4),vector(29,0),vector(24,-5),\
				 vector(20,0))',\
	lipso:hconcat(vector(0,-47),vector(7,-49),\
				  vector(16,-53),vector(7,-60),vector(0,-62))',\
	lipsi:hconcat(vector(7,-54),vector(0,-54)) ',\
	nose:hconcat(vector(0,-6),vector(3,-16),vector(6,-30),vector(0,-31))',\
	shape:hconcat(vector(0,44),vector(29,40),vector(51,22),\
				  vector(54,11),vector(52,-4),\
				  vector(46,-36),vector(38,-61),vector(25,-83),\
				  vector(0,-89))',\
	ear:hconcat(vector(60,-11),vector(57,-30))',\
	hair:hconcat(vector(72,12),vector(64,50),vector(36,74),\
				 vector(0,79))')
if (@oldvers) {
	# indices of features to be used in assignment
	@j_eye <- 1
	@j_iris <- 2
	@j_lipso <- 3
	@j_lipsi <- 4
	@j_nose <- 5
	@j_shape <- 6
	@j_ear <- 7
	@j_hair <- 8
}

@lipsiend <- 3 # index in lipso
@hairend <- 4  # index in shape
@earsta <- 5   # index in shape
@earend <- 6   # index in shape
@hair1 <- 1    # index in hair
@hair2 <- 2    # index in hair

@lipsoRefl <- run(4,1)  # lipso.refl.ind
@lipsiRefl <- 1         # lipsi.refl.ind
@noseRefl <- run(3,1)   # nose.refl.ind
@hairRefl <- run(3,1)   # hair.refl.ind
@shapeRefl <- run(8,1)  # shape.refl.ind<
@shapeXntN <- run(2,8)  # shape.xnotnull
@noseXntN <- run(2,3)   # nose.xnotnull

 #:5

 #6:
@result <- if (@keep){
	split(rep(0,@n)', compnames:delete(@xnames1,return:T))
} else {
	NULL
}

@ylim <- vector(-105,105)*1.2
@xlim <- 1.4*@ylim
for(@iface,1,@n){
 #7:
	@factors <- @xy[@iface,]
	@face <- @face_orig
 #:7

 #9:
	@lipso <- @face$lipso
	@lipsi <- @face$lipsi
	@m <- describe(@lipso[,2],mean:T)
	# adjust mouth height; default mouthsize=0.7
	@lipso[,2] <-\
		@m + (@lipso[,2] - @m)*(1 + @mods$mouthsize*@factors[4])
	@lipsi[,2] <-\
		@m + (@lipsi[,2] - @m)*(1 + @mods$mouthsize*@factors[4])
	# adjust mouth width
	@lipso[,1] <- @lipso[,1]*(1 + @mods$mouthsize*@factors[5])
	@lipsi[,1] <- @lipsi[,1]*(1 + @mods$mouthsize*@factors[5])
	# adjust curve of smile; default smile = 20
	@lipso[@lipsiend,2] <- @lipso[@lipsiend,2] + @mods$smile*@factors[6]

if (@oldvers) {
	@face[@j_lipso] <- delete(@lipso,return:T)
	@face[@j_lipsi] <- delete(@lipsi,return:T)
} else {
	@face$lipso <- delete(@lipso,return:T)
	@face$lipsi <- delete(@lipsi,return:T)
}
 #:9


 #10:
	@eye <- @face$eye
	@iris <- @face$iris
	@m <- describe(@eye[,2],mean:T)
	# adjust height of eyes; default eyesize = 0.7
	@eye[,2] <- @m + (@eye[,2] - @m)*(1 + @mods$eyesize*@factors[7])
	@iris[,2] <- @m + (@iris[,2] - @m)*(1 + @mods$eyesize*@factors[7])
	# adjust width of eyes
	@m <- describe(@eye[,1],mean:T)
	@eye[,1] <- @m + (@eye[,1] - @m)*(1 + @mods$eyesize*@factors[8])
	@iris[,1] <- @m + (@iris[,1] - @m)*(1 + @mods$eyesize*@factors[8])
if (@oldvers){
	@face[@j_eye] <- delete(@eye,return:T)
	@face[@j_iris] <- delete(@iris,return:T)
} else {
	@face$eye <- delete(@eye,return:T)
	@face$iris <- delete(@iris,return:T)
}
 #:10


 #11:
	@hair <- @face$hair
	@m <- min(@hair[,2])
	# adjust height of hair; default hairsize = 0.2
	@hair[,2] <- @m + (@hair[,2] - @m)*(1 + @mods$hairsize*@factors[9])
	@m <- 0
	# adjust width of hair
	@hair[,1] <- @m + (@hair[,1] - @m)*(1 + @mods$hairsize*@factors[10])
	@m <- 0
	# adjust styling of hair; default hairstyle = 50
	@hair[vector(@hair1,@hair2),2] <- \
		@hair[vector(@hair1,@hair2),2] + @mods$hairstyle*@factors[11]
if (@oldvers){
	@face[@j_hair] <- delete(@hair,return:T)
} else {
	@face[@j_hair] <- delete(@hair,return:T)
}
 #:11


 #12:
	@nose <- @face$nose
	@m <- describe(@nose[,2],mean:T)
	# adjust height of nose; default noseheight = 0.7
	@nose[,2] <- @m + (@nose[,2] - @m)*(1 + @mods$noseheight*@factors[12])
	# adjust width of nose; default nosewidth = 1.0
	@nose[@noseXntN,1] <-\
		@nose[@noseXntN,1]*(1 + @mods$nosewidth*@factors[13])
if (@oldvers) {
	@face[@j_nose] <- delete(@nose,return:T)
} else {
	@face$nose <- delete(@nose,return:T)
}
 #:12


 #13:
	@ear <- @face$ear
	@m <- describe(@face$shape[vector(@earsta,@earend),1],mean:T)
	# adjust width of ears; default earsize = 0.7
	@ear[,1] <- @m + (@ear[,1] - @m)*(1 + @mods$earsize*@factors[14])
	@m <- min(@ear[,2])
	# adjust height of ears
	@ear[,2] <- @m + (@ear[,2] - @m)*(1 + @mods$earsize*@factors[15])
if (@oldvers) {
	@face[@j_ear] <- delete(@ear,return:T)
} else {
	@face$ear <- delete(@ear,return:T)
}

 #:13


 #8:
	for(@i,1,ncomps(@face)) {
		@comp <- @face[@i]
		# adjust height of face; default facesize = 0.2
		@comp[,2] <- @comp[,2]*(1 + @mods$facesize*@factors[1])
		# adjust width of face
		@comp[,1] <- @comp[,1]*(1 + @mods$facesize*@factors[2])
		@J1 <- @comp[,1] <= 0
		if (sum(@J1) > 0) {
			@comp[@J1,1] <- 0
		} elseif (sum(!@J1) > 0) {
			@J2 <- @comp[,2] <= -30
			if (sum(@J2) > 0) {
			# adjust shape of face; default faceshape = 0.2
				@tmp <- @comp[@J2,1] +\
					(@comp[@J2,2] + 50)*@mods$faceshape*\
						sin(-1.5*@factors[3],radians:T)
				@tmp[@tmp < 0] <- 0
				@comp[@J2,1] <- delete(@tmp,return:T)
			}
		}
		@face[@i] <- @comp
	}
	delete(@comp,@i,@J1)

 #:8

 #14:
	if (!ismacro(@invert)) {
		@invert <- macro("(\$1) * vector(-1,1)'")
	}
	# Build complete face
	@face_obj <- structure(\
		eyer:@face$eye, \
		eyel:@invert(@face$eye), \
		irisr:@face$iris, \
		irisl:@invert(@face$iris), \
		lipso:vconcat(@face$lipso,@invert(@face$lipso[@lipsoRefl,])), \
		lipsi:vconcat(@face$lipso[@lipsiend,],@face$lipsi,\
					 @invert(@face$lipsi[@lipsiRefl,]),\
					 @invert(@face$lipso[@lipsiend,])), \
		earr:vconcat(@face$shape[@earsta,],@face$ear,@face$shape[@earend,]), \
		earl:@invert(vconcat(@face$shape[@earsta,],@face$ear,@face$shape[@earend,])), \
		nose:vconcat(@face$nose,@invert(@face$nose[@noseRefl,])), \
		hair:vconcat(@face$shape[@hairend,],@face$hair,@invert(@face$hair[@hairRefl,]),\
					@invert(@face$shape[@hairend,])), \
		shape:vconcat(@face$shape,@invert(@face$shape[@shapeRefl,])))
 #:14

 #15:
	if (!ismacro(interpspline)) {
		getmacros(interpspline,silent:T)
	}
	@open <- vector(6,7,8,10)
	@closed <- run(ncomps(@face_obj))[-@open]
	@Xcoords <- @Ycoords <- split(rep(0,ncomps(@face_obj))')
	@k <- 0
	for(@jcomp,@open) {
		@x <- @face_obj[@jcomp][,1]
		@y <- @face_obj[@jcomp][,2]
		@Xcoords[@jcomp] <- \
			vector(interpspline(run(length(@x)),@x,npoints:@npts)[,2],?)
		@Ycoords[@jcomp] <- \
			vector(interpspline(run(length(@y)),@y,npoints:@npts)[,2],?)
	}
	for(@jcomp,@closed) { # will use periodic spline
		@x <- @face_obj[@jcomp][,1]
		@y <- @face_obj[@jcomp][,2]
		@Xcoords[@jcomp] <- \
			vector(interpspline(run(length(@x)),@x,npoints:@npts)[,2],?)
		@Ycoords[@jcomp] <- \
			vector(interpspline(run(length(@y)),@y,npoints:@npts)[,2],?)
	}
	@Xcoords <- vector(@Xcoords)
	@Ycoords <- vector(@Ycoords)
	#print(xcoords:paste(@Xcoords,missing:"?"))
	#print(ycoords:paste(@Ycoords,missing:"?"))
	if (@keep) {
		# add two isolated corner points so all faces have same scale
		@Xcoords <- vector(@Xcoords,@xlim[1],?,@xlim[2])
		@Ycoords <- vector(@Ycoords,@ylim[1],?,@ylim[2])
		@result[@iface] <- hconcat(@Xcoords,@Ycoords)
	}
	
	if (@draw) {
		if (!@panels) {
			lineplot(@Xcoords,@Ycoords,title:@xnames[@iface],\
				   xlab:paste("Face derived from case", @iface),\
				   $K,xaxis:F,yaxis:F,xticks:NULL,yticks:NULL,\
				   xmin:@xlim[1],xmax:@xlim[2],ymin:@ylim[1],ymax:@ylim[2])
		} else {
			if (!ismacro(panelplot)) {
				getmacros(panelplot,silent:T)
			}
			panelplot(@Xcoords,@Ycoords,pos:@iface,\
					hstrips:@nr,vstrips:@nc,lines:T,\
					xaxis:F,yaxis:F,\
					xmin:@xlim[1],xmax:@xlim[2],ymin:@ylim[1],ymax:@ylim[2],\
					show:F,gridlines:@iface==1 && @gridlines,add:@iface>1)
			if (@xnames[@iface] != "") {
				@i <- ceiling(@iface/@nc)
				@j <- 1 + (@iface - 1) %% @nc
				addstrings(@j - 1 + @horiz,@nr - @i + @vert, \
						   string:@xnames[@iface], just:@just, show:F)
			}
		}
	}
}
if (@oldvers) {
	delete(@j_eye,@j_iris,@j_lipso,@j_lipsi,@j_nose,@j_shape,@j_ear,@j_hair)
}

if (@panels) {
	showplot(notakey:NULL,$K,xlab:"Cases in order left to right by rows",\
			 title:"Faces plot of $1")
	delete(@vert,@horiz,@just)
}
 #:15
delete(@iface,@ylim,@xlim,@Xcoords,@Ycoords,@jcomp,@x,@y,@npts,@panels,\
	@xnames,@gridlines,@xy,@mods,@mm,@n,@nr,@nc,\
	@face_orig,@face,@face_obj,@n_char,\
	@lipsiend,@hairend,@earsta,@earend,@hair1,@hair2,\
	@lipsoRefl,@lipsiRefl,@noseRefl,@hairRefl,@shapeRefl,@shapeXntN,\
	@noseXntN,@m,@factors)

delete(@result,return:T)
%faces%

===> chiqqplot <===
chiqqplot       MACRO DOLLARS
) Macro to make a chisquare or sqrt(chisquare) QQ plot of generalized
) distance of residuals from 0
) usage:
)   chiqqplot(dsq, df [,sqrt:T] [, symbol:sym])
) or
)   chiqqplot(y [,df] [,quadform:Q] [,center:cent] [,sqrt:T] \
)                [, symbol:sym] [graphics keyword phrases)
)     dsq               REAL vector or matrix of squared distances >= 0
)     df                REAL scalar > 0
)     y                 REAL n by p data matrix with no MISSING values
)                       from which a vector dsq of squared generalized
)                       distances from is computed and plotted
)     Q                 REAL positive definite p by p matrix, default
)                       is S = sample covariance matrix of y
)     cent              REAL length p row or column vector with no
)                       MISSING values; default is ybar = sample
)                       mean vector
)     sqrt:T            plot square roots
)     sym               integer or CHARACTER vector or matrix
) In the 2nd usage
)   dsq[i] = (y[i,] - cent') %*% solve(Q) %*% (y[i,] - cent')'
)   The default is Q = S = the sample variance matrix.
)   The default value for df is ncols(y)
)
) With df but without quadform:Q, argument 1 is treated as a vector
) of distances or matrix whose columns are considered as vectors of
) distances, all to be plotted on the same grapn.
)
) With quadform:Q or without df, argument 1 is treated as a data matrix
) or vector.
)
) Ordered values of dsq or each column of dsq are plotted against
) equally spaced quantiles of chi-squared on df degrees of freedom.
)
) With 'sqrt:T', sqrt(dsq) is plotted against sqrt(chi-squared quantiles)
)
) With symbol:0, when ncols(d) = 1, the symbols are the case numbers, and
) when ncols(d) > 1, the symbols are the column numbers.  Otherwise,
) symbol:sym works as for plot().
)
) You can use the usual graphics keywords such as 'xlab', 'ylab' and
) 'title'
)
))Written 000601 by C. Bingham, kb@umn,edu
))040102 Added keywords quadform, center
#$S(d, df [, sqrt:T] [,symbols:T] [,graphics keyword phrases])
#$S(y [,df] [,quadform:Q] [,center:cent] [,sqrt:T] \
#   [, symbol:sym] [graphics keyword phrases)
if ($v > 2){
	error("usage: $S(y [,df] [keyword phrases]) or $S(d, df [keyword phrases])")
}
@y <- argvalue($1, "argument 1","nonmissing real matrix")
@df <- if ($v > 1) {
	argvalue($2,"df","positive number")
} else {
	NULL
}
@Q <- keyvalue($K,"quadf*","nonmissing real square")
@cent <- keyvalue($K,"cent*","real matrix")
if (!isnull(@df) && isnull(@Q) && isnull(@cent)){
	rename(@y,@dsq)
	if (min(vector(@dsq)) < 0) {
		error("at least 1 distance < 0")
	}
} else { # arg 1 is data matrix
	@p <- ncols(@y)
	if (isnull(@df)) {
		@df <- @p
	}
	if (isnull(@cent)) {
		@cent <- describe(@y,mean:T)
	} else {
		if (nrows(@cent) == 1) {
			@cent <- @cent'
		} elseif (!isvector(@cent)) {
			error("value of 'center' not row or column vector")
		}
		if (length(@cent) != @p) {
			error("length(center) != ncols(y)")
		}
	}
	
	if (!isnull(@Q)) {
		if (max(abs(vector(@Q - @Q'))) >= 1e-7) {
			error("quadratic form matrix not symmetric")
		}
		if (nrows(@Q) != @p) {
			error(paste("quadratic form matrix not",@p,"by",@p))
		}
	} else {
		@Q <- tabs(@y, covar:T)
	}
	@vals <- eigenvals(@Q)
	if (anytrue(@vals[1] <= 0, @vals[@p]/@vals[1] <= 1e-9)) {
		error("variance or quadratic form matrix not positive definite")
	}
	delete(@vals)
	@y <-- @cent'
	@dsq <- vector((@y * (@y %/% @Q)) %*% rep (1,@p))
	delete(@y, @p)
}
delete(@cent,@Q)

@n <- nrows(@dsq)
@m <- ncols(@dsq)

@sqrt <- keyvalue($K,"sqrt","TF",default:F)
@ch <- vector("\1","\2","\3","\4","\5","\6","\7","\10")[1 + run(0,@m-1)%%8]
@symbols <- keyvalue($K,"symbol*","matrix",default:delete(@ch,return:T))

if (alltrue(!ischar(@symbols),!isreal(@symbols, integer:T, nonneg:T))) {
	error("Value of 'symbols' not CHARACTER or non-negative integer")
}
if (isreal(@symbols)) {
	if (max(@symbols) > 999) {
		error("largest symbol number >= 1000")
	}
	if (issalar(@symbols) && @symbols[1] == 0) {
		@symbols <- if (@m == 1){
			grade(@dsq)
		}else{
			run(@m)'
		}
	}
}
@dsq <- sort(@dsq)
@vals <- invchi((run(@n) - .5)/@n,@df)
if (@sqrt){
	@dsq <- sqrt(@dsq)
	@vals <- sqrt(@vals)
	@what <- "Chi"
	@ylab <- "Distances"
}else{
	@what <- "Chi-squared"
	@ylab <- "Squared distances"
}
@xlab <- paste(@what,"(",@df,") probability points",sep:"")
@title <- paste(@what, "Q-Q plot of $1")

plot(delete(@vals,return:T), delete(@dsq,return:T),\
	symbols:delete(@symbols,return:T), $K,\
	xlab:delete(@xlab,return:T),ylab:delete(@ylab,return:T),\
	title:delete(@title,return:T),xmin:0, ymin:0)
%chiqqplot%

===> cluscritplot <===
cluscritplot  macro  dollars
) Macro to make plots to aid in choosing the number of clusters from
) the vector of criterion values computed by cluster().
)
) cluscritplot() plots against stage values of and/or changes in
) the criterion (inter-cluster distance of clusters merged at that stage)
)
) Usages:
)  cluscritplot(Criterion [, N] [,changes:T] [,values:T] [,quiet:T]\
)               [, graphics keywords])
)  cluscritplot(str [,N] [,changes:T] [,values:T] [,quiet:T]\
)               [, graphics keywords])
)   Criterion    positive vector, length > 1, either monotonically
)                decreasing or increasing
)   N            positive integer = sample size; if omitted, 
)                N = length(Criterion) + 1 is used
)   str          structure with component 'criterion' and component
)                'classes' like the value of cluster(...,keep:"all") or
)                cluster(...,keep:vector("criterion","classes"))
)                Criterion = str$criterion and N = nrows(str$classes)
)                unless N is provided as argument 2
)   changes:T    plot changes in criterion (default without values:T);
)                can be used with values:T
)   values:T     plot the criterion itself (default with changes:F);
)                can be used with changes:T
)   quiet:T      suppress any warning messages
)
) Criterion is reordered, if necessary, in increasing order.
)
) With changes:T, it plots against stage number the first differences
) Criterion[2]-Criterion[1], ..., Criterion[m]-Criterion[m-1], where
) m = length(Criterion).  Points are labelled with the number of
) clusters before the change, the last point being labeled with "2".
)
) With values:T, Criterion is plotted against stage number with final
) stage numbered n-1 Points are labeled with the number of clusters
) before the merge at that stage merge.
)
))031205 written by C. Bingham, kb@umn.edu
))040102 Minor mainly cosmetic modifications
))040912 Correct minor bug
#$S(Criterion [, N] [,changes:T] [,values:T] [,quiet:T] [, graphics keywords])
#$S(str, [, N] [,changes:T] [,values:T] [,quiet:T] [, graphics keywords])
@arg1 <- argvalue($1,"criterion or structure")
@n <- NULL

@crit <- if (isstruc(@arg1)) {
	if (match("criterion",compnames(@arg1),0) == 0) {
		error("No 'criterion' component in structure $1")
	}
	if (match("classes",compnames(@arg1),0) != 0) {
		@n <- nrows(@arg1$classes)
	}
	argvalue(@arg1$criterion,"component criterion","positive vector")
} else {
	argvalue(@arg1,"criterion","positive vector")
}
delete(@arg1)

@m <- length(@crit)
if (@m < 2) {
	error("criterion vector has only 1 element")
}

@quiet <- keyvalue($K,"quiet","TF",default:F)

@changes <- keyvalue($K,"change*","TF")
@values <- keyvalue($K,"value*","TF")
if (isnull(@changes) && isnull(@values)) {
	@changes <- T
	@values <- F
} elseif (isnull(@changes)) {
	@changes <- !@values
} elseif (isnull(@values)) {
	@values <- !@changes
} elseif (!@changes && !@values) {
	error("can't have both changes:F and values:F")
}

@jumps <- movavg(1,@crit)[-1] # first differences
if (min(@jumps) < 0) {
	if (max(@jumps) > 0) {
		error("Criterion vector not monotonic")
	}
	@crit <- reverse(@crit)
	@jumps <- -reverse(@jumps)
}
if ($v > 1) {
	@n <- argvalue($2,"n","positive count")
} elseif (isnull(@n)) {
	if (!@quiet) {
		print(paste("WARNING: no sample size available; using m + 1 =",@m+1))
	}			
	@n <- @m + 1
}

if (@values) {
	plot(@n-@m,@crit,$K,symbols:run(@m,1),\
	     xlab:"Stage",ylab:"Criterion",ymin:0,ymax:1.05*max(@crit),\
	     lines:T, title:"Merging criterion vs stage")
}
if (@changes) {
	plot(@n-@m+1,@jumps,$K,symbols:run(@m,2),\
	     xlab:"Stage", ylab:"Change",ymin:0, ymax:1.05*max(@jumps),\
	     lines:T, title:"Change in criterion from preceding stage vs stage")
}
delete(@changes,@values,@quiet,@n,@m,@crit,@jumps)
%cluscritplot%

_E_N_D_O_F_M_A_C_R_O_S_ Marker for end of macros and data

Help file for mvplots.mac for MacAnova
(C) 2003 by Gary W. Oehlert and Christopher Bingham
Updated 040912 CB

!!!! Starting marker for message of the day
!!!! Ending marker for message of the day

???? Starting marker for list of up to 32 comma/newline separated keys
Cluster Analysis
Multivariate
Plotting
Splines
???? Ending marker for keys

====andrewsplot()#multivariate, plotting
%%%%
andrewsplot(x [,rescale:T or stand:T] [,npoints:np] [,graph keys])
  REAL data matrix x with no MISSING values, integer np > 1
andrewsplot(x, group [,linetypes:ltype] [,rescale:T or stand:T]\
  [,npoints:np] [,graph keywords]), positive integer vector group
%%%%
@@@@usage
andrewsplot(x) draws an Andrews plot of the multivariate data in a n
by p REAL matrix x with no MISSING data.  Specifically, for each row
x[i,], a trigonometric function f(t) is computed and plotted as a line
graph against equally spaced values of t between -.525 and .525.

The function for case i is

  f(t) = x[i,1]/sqrt+x[i,2]*cos(t)+x[i,3]*sin(t)+x[i,4]*cos(2*t) ...

where the arguments to cos() and sin() are in cycles.  The last term is
x[i,p]*sin(m*t) or x[i,p]*cos(m*t) depending on whether p is even or
odd, where m = floor(p/2)

andrewsplot(x, group [,linetypes:ltypes]), where group is a vector of
positive integers with length(group) = nrows(x), does the same except a
different line type is used for each value of group.  When
linetypes:ltypes is an argument, ltypes must be a vector of positive
integers with length(ltypes) >= max(group).

@@@@keywords
With either usage, you can include 'npoints:np' and/or either of
'stand:T' or 'rescale:T' as additional arguments.

With 'npoints:np', where np is an integer > 2, the functions are
computed and plotted at np equally spaced values of t.  The default is
np = 35.

With 'rescale:T', x is rescaled as (x - min(x))/(max(x) - min(x)) before
f(t) is computed.

With 'stand:T', x is standardized as (x - xbar)/stddev before f(t) is
computed.

@@@@graphics_keywords
As additional arguments, you can use graphics keyword phrases add:T,
'show:F', and 'window:n' plus any graphics keyword phrases that can be
used with showplot().  See 'graph_keys', showplot().

@@@@see_also
See also parcoordplot(), lineplot(), rowplot().

====cluscritplot()#multivariate, plotting
%%%%
cluscritplot(Criterion [,N] [,changes:T or values:T] [,quiet:T]\
  [,graphics keywords]), positive increasing or decreasing REAL vector
  Criterion, N > 0 integer
cluscritplot(str [,N] [,changes:T or values:T] [,quiet:T]\
  [,graphics keywords]), str = structure(criterion:Criterion,
  [,classes:Classes]), Classes vector of positive integers
%%%%
@@@@introduction#Introduction
cluscritplot() makes graphs that help in choosing the number of clusters from
the vector of criterion values computed by cluster().  The criterion values
are the distances between the two clusters that are merged at each stage.

cluscritplot() plots against stage number criterion values and/or
stage-to-stage changes in the criterion

@@@@usage#Usage
There are two essentially equivalent usages, one in which the criterion
values are provided in a vector and the other in which they are provided
as a component of a structure as is returned by cluster() with argument
'keep:vector("criterion","classes")'.

cluscritplot(Crit, N) and cluscritplot(Crit, N, changes:T) plots the
changes d[i] = Crit[i+1] - Crit[i] against i + N - M, i = 1,...,
M-1, where M = length(Crit).  Crit must be an increasing (d[i] >= 0)
vector and N > 0 an integer = size of sample being cluster.  The
Final point is at N-1, the number of the stage just before the merge of
the final two clusters.

Crit can also be in decreasing order, in which case reverse(Crit) is
used.  If N is omitted, d[i] is plotted is plotted against i = 2, ...,
M.

The graph is a line connected plot, with integer plotting symbols with
the final point labeled "2".  When there is a sharp peak at a point
labeled with integer k, it suggests that k might be a candidate for the
number of clusters.

cluscritplot(Str) and cluscritplot(Str, changes:T), where Str =
structure(classes:Classes,criterion:Crit) as is returned by
cluster(y,..., keep=vector("criterion","classes")) is equivalent to
cluscritplot(Str$criterion,nrows($Str$classes).

cluscritplot(Crit, N, values:T) and cluscritplot(Str, values:T) make
similar plots except Crit[i] is plotted against i + N - M - 1, i =
1,...,M (1, ..., M when N is omitted).  The points are labeled so that
the final point is labeled 1.  When there is a sharp rise in the plot at
the point labelled k, it suggests that k might be a candidate for the
number of clusters.

The usual graphics keywords such as 'title' and 'linetype' may be used.

@@@@keyword_quiet#Keyword 'quiet'
When keyword phrase 'quiet:T' is an argument, certain warning messages
are suppressed.

@@@@see_also#References
See also cluster(), 'graph_keys', 'graphs'.

====faces()#multivariate, plotting
%%%%
faces(Y [,whichvars:varnos] [,byrow:T] [,nrows:nr] [,ncols:nc]  \
  [,min:Min] [,max:Max] [,scale:F] [,panel:F] [,keep:T or draw:T] \
  [,labels:facelabs] [,labelplace:lp] [,gridlines:T] [,npoints:npts] \
  [,mods:Mods] [graph keywords]), Y REAL nonMISSING matrix, nr > 0, nc >
  0 integers, Mods a structure with positive scalar components,
  CHARACTER vector facelabs, CHARACTER scalar lp, integer npts > 0
%%%%
@@@@introduction#Introduction
faces() is a macro to draw Chernoff type faces from multivariate data.

@@@@usage#Usage
faces(Y) draws N Chernoff faces, one for each row in N by p data matrix
Y.  Each face is framed in its own pane.  You can use most of the
standard graphics keywords.

faces(Y, nr:nrows, nc:ncols) does the same, except the faces are
arranged in nrows rows of up to ncols faces each.  It is an error if
nrows*ncols < N.

faces(Y, whichvars:varnos, ...) is the same as faces(Y[,varnos], ...),
where varnos is a vector of integers between 1 and p = ncols(Y).

faces(Y, byrow:T ...) is the same as faces(Y' ...), that is columns of Y
represent cases and rows represent variables.

By default, before computing the faces, Y is rescaled to
X = 2*(Y - Min)/(Max - Min) - 1, Min = min(Y), Max = max(Y), so that all
values are between -1 and 1.  You can suppress this by including
'scale:F' as an argument.  This will usually be a disaster unless Y has
previously been suitably scaled.

Occasionally it may be desired to modify the scaling by providing
vectors Min and/or Max by keyword phrases min:Min and max:Max where Min
and Max are row or column vectors.  This option is important if you are
plotting faces of different sets of cases in different windows.  In
that case you should use the same Min and Max in every window to ensure
that the faces in different windows are scaled the same.

@@@@keep_and_draw keywords#Keywords 'keep' and 'draw'
facecoords <- faces(Y, keep:T) does not draw faces.  Instead it sets
facecoords to structure(face_1, face_2, ...), where face_j is a matrix
with two columns containing x and y coordinates for the face computed
from Y[j,].

facecoords <- faces(Y, keep:T, draw:T), does the same, except faces are
also drawn.

@@@@gridlines_keyword#Keyword 'gridlines'
faces(Y, gridlines:F ...) draws faces without lines between panes.

@@@@panels_keyword#Keyword 'panels'
faces(Y, panels:F) draws each face in a separate window.  It is an error
if there are not at least nrows(Y) empty graphics windows.

@@@@face_labels#Labelling faces
By default, faces are labeled 1, 2, ..., nrows(Y).  If Facelabs is a
CHARACTER vector of length N, faces(Y,labels:Facelabs) labels the face
associated with Y[i,] with Facelabs[i].

You can control the position of the labels using keyword phrase
'labelplace:labpos', where labpos = "vh" where v is 'b' (bottom), 'm'
(middle) or 't' (top) specifying the vertical position and h is 'l'
(left), 'c' (center) or 'r' (right) specifying the horizontal position.
labelplace:"mc" is not a good idea.

@@@@npoints_keyword#Keyword 'npoints'
Each feature is drawn with an interpolating spline computed at M points,
where 10 <= M <= 40, decreasing as N increases (M = 40 with panels:F).
You can overwrite this with keyword phrase 'npoints:M', where M > 0 is
an integer.

@@@@mods_keyword#Keyword 'mods'
If you don't like the mapping of the scaled Y values into features, you
can modify it by including mods:Mods as an argument where Mods is a
structure with positive scalar components with names chosen from
'facesize' [.2], 'faceshape' [.2], 'mouthsize' [.7], 'smile' [20],
'eyesize' [.7], 'hairsize' [.2], 'hairstyle' [50], 'noseheight' [.7],
'nosewidth' [1], or 'earsize' [.7] used to replace default values (in
[...]) used to translate variables into feature characteristics; see
below

@@@@shaping_features#Shaping features
There are 15 features modified by the data.  In the list below, x and y
to the left of '=' are the plotted coordinates; x and y to the right of
'=' are the coordinates of the base features.  X[1], ..., X[15] are the
(usually) scaled values of Y[i,1], Y[i,2], ..., Y[i,p], p = ncols(Y).
When p < 15, the variables are repeated cyclically.  If p > 15, extra
variables are ignored

         Features              Mapping         
     height of face   y = y*(1 + facesize*X[1]), all features
     width of face    x = x*(1 + facesize*X[2]), all features
     shape of face    x = x + (y + 50)*faceshape*sin(-1.5*X[3])
     height of mouth  y = m + (y-m)*(1 + mouthsize*X[4]), m = mean
     width of mouth   x = x*(1 + mouthsize*X[5])
     curve of smile   x = x + smile*X[6])
     height of eyes   y = m + (y-m)*(1 + eyesize*X[7]), m = mean
     width of eyes    x = m + (x-m)*(1 + eyesize*X[8]), m = mean
     height of hair   y = m + (y-m)*(1 + hairsize*X[9]), m = min
     width of hair    x = m + (x-m)*(1 + hairsize*X[10]), m = 0
     styling of hair  x = x + hairstyle*X[11])
     height of nose   y = m + (y-m)*(1 + noseheight*X[12]), m = mean
     width of nose    x = x*(1 + nosewidth*X[13])
     width of ears    x = m + (x-m)*(1 + earsize*X[14]), m = mean
     height of ears   y = m + (y-m)*(1 + earsize*X[15]), m = min

@@@@references#References
faces() was translated and modified from R code of Dr. Hans Peter Wolf
(pwolf@wiwi.uni-bielefeld.de), downloaded from
http://www.wiwi.uni-bielefeld.de/~wolf/software/R-wtools/faces/.  It
computes a spline interpolant with macro interpspline(), which is also
adapted from R code of Hans Peter Wolf, based on Ralston and Rabinowtz
(1985), A First Course in Numerical Analysis.

@@@@see_also
See also andrewsplot(), parcoordplot(), interpspline().

====interpspline()#splines
%%%%
result <- interpspline(a, y [, npoints:m] [,plot:T]), nonMISSING REAL
  vectors, integer m > 0
%%%%
@@@@introduction#Introduction
Macro to compute spline interpolant.  It is adapted from R code of Hans
Peter Wolf, based on Ralston and Rabinowitz (1985), A First Course in
Numerical Analysis, pp. 76ff.

@@@@usage#Usage
Usage information is taken from comments assocviated with interpspline()
and is preliminary.
Usage
 result <- interpspline(a, y)
  a         nonMISSING REAL vector with a[i] > a[i-1]
  y         nonMISSING REAL vector with n = length(y) = length(a)
  plot:T    make line plot
  result    hconcat(x, Sx), x = a[1] + (a[n]-a[1])*run(0,m-1)/(m-1),
            Sx = REAL vector of length m containing spline values
            at equally spaced values x.

A cubic spline is determined to interpolate the points (a[1],y[1]), ...,
(a[n],y[n]).  Then the values of the spline are computed at m equally
spaced points x[1], x[2], ..., x[n], with x[1] = a[a] and x[m] = a[n].
The default value for m is 50.

Thus, say, lineplot(result[,1], result[,2]) will plot a cubioc spline
curve that exactly goes through the points (a[1],y[1]), ..., (a[n],y[n]).

@@@@keywords#Keywords
result <- interpspline(a, y, npoints:M), M > 0 an integer does the same
but computes the spline at m = M points.

result <- interpspline(a, y, plot:T [, npoints:M] [graphics keywords])
does the same, except the spline curve is drawn with superimposed
symbols at the points being interpolated.

====parcoordplot()#multivariate, plotting
%%%%
parcoordplot(x [,rescale:F or stand:T] [,graph keys]), REAL matrix x
parcoordplot(x, group [,linetypes:ltype] [,rescale:F or stand:T]\
   [,graph keys]), group and ltype vectors of positive integers
%%%%
parcoordplot(x) draws a parallel coordinates plot of rescaled data from
a n by p REAL matrix x.  The rescaled version of x is y =
(x - min(x)/(max(x) - min(x)).  Specificaly, for each case, a line plot
is drawn of y[i,] vs i, i = 1, ...,p.

parcoordplot(x, group [,linetypes:ltypes] [,npoints:np]), where group is
a vector of positive integers with length(group) = nrows(x), does the
same except a different line type is used for each value of group.  When
linetypes:ltypes is an argument, ltypes must be a vector of positive
integers with length(ltypes) >= max(group).

@@@@keywords
With either usage, you can include either of 'stand:T' or 'rescale:F' as
additional arguments.

With 'stand:T', x is standardized to y = (x - xbar)/stddev instead of
being rescaled.

With 'rescale:F' and without 'stand:T', x is neither rescaled or
standardized.

@@@@graphics_keywords
As additional arguments, you can use graphics keyword phrases add:T,
'show:F', and 'window:n' plus any graphics keyword phrases that can be
used with showplot().  See 'graph_keys', showplot().

@@@@see_also
See also andrewsplot(), faces(), lineplot(), rowplot().


_E_O_F_#This should be the last line, an internal End Of File marker