Skip to content

Commit bb3c7f2

Browse files
friendlyfriendly
committed
test plotEqn
1 parent 52699f1 commit bb3c7f2

File tree

5 files changed

+347
-2
lines changed

5 files changed

+347
-2
lines changed

‎NEWS.md

+1
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@
44
- Add `papers/matlib-useR-2016.pdf` to avoid bad URL
55
- Consolidate options for `print.latexMatrix`
66
- Fix bug in `print.latexMatrix(sparse=TRUE)`
7+
- `plotEqn()` gains a `...` to pass other graphical parameters
78

89
# matlib 1.0.0
910

‎dev/plotEqn-test.R

+22
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,22 @@
1+
# show dual relation between points & lines
2+
3+
library(matlib)
4+
5+
A <- matrix(c( 1, 2, 0,
6+
-1, 2, 1), 3, 2) |>
7+
print()
8+
9+
10+
b <- c(2, 1, 1)
11+
12+
showEqn(A, b, vars = c("x", "y"), simplify = TRUE)
13+
14+
plotEqn(A, b, vars = c("x", "y"),
15+
cex.lab = 2,
16+
solution = list(pch = 16))
17+
18+
# try to change the labels: doesn't work
19+
plotEqn(A, b, vars = c("x", "y"),
20+
labels = c("y = x - 2",
21+
"y = 1/2 - x",
22+
"y = 1"))

‎dev/plotEqn.R

+319
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,319 @@
1+
#' Plot Linear Equations
2+
#'
3+
#' Shows what matrices \eqn{A, b} look like as the system of linear equations, \eqn{A x = b} with two unknowns,
4+
#' x1, x2, by plotting a line for each equation.
5+
#'
6+
#' @param A either the matrix of coefficients of a system of linear equations, or the matrix \code{cbind(A,b)}.
7+
#' The \code{A} matrix must have two columns.
8+
#' @param b if supplied, the vector of constants on the right hand side of the equations, of length matching
9+
#' the number of rows of \code{A}.
10+
#' @param vars a numeric or character vector of names of the variables.
11+
#' If supplied, the length must be equal to the number of unknowns in the equations, i.e., 2.
12+
#' The default is \code{c(expression(x[1]), expression(x[2]))}.
13+
#' @param xlim horizontal axis limits for the first variable
14+
#' @param ylim vertical axis limits for the second variable; if missing, \code{ylim} is calculated from the
15+
#' range of the set of equations over the \code{xlim}.
16+
#' @param col scalar or vector of colors for the lines, recycled as necessary
17+
#' @param lwd scalar or vector of line widths for the lines, recycled as necessary
18+
#' @param lty scalar or vector of line types for the lines, recycled as necessary
19+
#' @param axes logical; draw horizontal and vertical axes through (0,0)?
20+
#' @param labels logical, or a vector of character labels for the equations; if \code{TRUE}, each equation is labeled
21+
#' using the character string resulting from \code{\link{showEqn}}, modified so that the
22+
#' \code{x}s are properly subscripted.
23+
#' @param ... Other arguments passed to \code{plot}
24+
#' @param solution logical: should the solution points for pairs of equations be marked? This can also be a list
25+
#' giving graphical parameters for the solution points.
26+
#' @return nothing; used for the side effect of making a plot
27+
#'
28+
#' @author Michael Friendly
29+
#' @references Fox, J. and Friendly, M. (2016). "Visualizing Simultaneous Linear Equations, Geometric Vectors, and
30+
#' Least-Squares Regression with the matlib Package for R". \emph{useR Conference}, Stanford, CA, June 27 - June 30, 2016.
31+
#' @importFrom graphics abline lines plot text points
32+
#' @export
33+
#' @seealso \code{\link{showEqn}}, \code{vignette("linear-equations", package="matlib")}
34+
35+
#' @examples
36+
#' # consistent equations
37+
#' A<- matrix(c(1,2,3, -1, 2, 1),3,2)
38+
#' b <- c(2,1,3)
39+
#' showEqn(A, b)
40+
#' plotEqn(A,b)
41+
#'
42+
#' # inconsistent equations
43+
#' b <- c(2,1,6)
44+
#' showEqn(A, b)
45+
#' plotEqn(A,b)
46+
47+
plotEqn <- function(A, b, vars, xlim, ylim,
48+
col=1:nrow(A),
49+
lwd=2, lty=1,
50+
axes=TRUE, labels=TRUE,
51+
solution=TRUE,
52+
...
53+
) {
54+
55+
if (!is.numeric(A) || !is.matrix(A)) stop("A must be a numeric matrix")
56+
if (missing(b)) {
57+
b <- A[ , ncol(A)] # assume last column of Ab
58+
A <- A[ , -ncol(A), drop=FALSE] # remove b from A
59+
}
60+
if (ncol(A) != 2) stop("plotEqn only handles two-variable equations. Use plotEqn3d for three-variable equations.")
61+
62+
if (missing(vars)) vars <- c(expression(x[1]), expression(x[2])) # paste0("x", 1:ncol(A))
63+
64+
neq <- nrow(A)
65+
66+
# establish x-axis limits and preliminary y-axis limits based on equation intersections
67+
68+
if (missing(xlim) || missing(ylim)) {
69+
if (neq == 1){
70+
if (missing(xlim)) xlim <- c(-4, 4)
71+
ylim.0 <- NULL
72+
intersections <- NULL
73+
} else {
74+
intersections <- matrix(NA, nrow=neq*(neq - 1)/2, ncol=2)
75+
colnames(intersections) <- c("x", "y")
76+
k <- 0
77+
for (i in 1:(neq - 1)) {
78+
for (j in (i + 1):neq) {
79+
k <- k + 1
80+
x <- try(solve(A[c(i, j), ], b[c(i, j)]), silent=TRUE)
81+
if (!inherits(x, "try-error")) intersections[k, ] <- x
82+
}
83+
}
84+
if (missing(xlim)) {
85+
xlim.0 <- if (length(unique(signif(intersections[, 1]))) != 1){
86+
c(-1, 1) + range(intersections[ , 1], na.rm=TRUE)
87+
} else c(-5, 5) + intersections[1, 1]
88+
xlim <- if (!any(is.na(xlim.0))) xlim.0 else c(-4, 4)
89+
}
90+
if (missing(ylim)) {
91+
ylim.0 <- if (length(unique(signif(intersections[, 2]))) != 1){
92+
c(-1, 1) + range(intersections[ , 2], na.rm=TRUE)
93+
} else c(-5, 5) + intersections[1, 2]
94+
if (any(is.na(ylim.0))) ylim.0 <- NULL
95+
}
96+
}
97+
}
98+
99+
# set values for horizontal variable
100+
x <- seq(xlim[1], xlim[2], length=10)
101+
102+
if (length(col) < neq) col <- rep_len(col, length.out=neq)
103+
if (length(lwd) < neq) lwd <- rep_len(lwd, length.out=neq)
104+
if (length(lty) < neq) lty <- rep_len(lty, length.out=neq)
105+
106+
if (missing(ylim)) {
107+
ylim <- ylim.0
108+
for (i in 1:neq) {
109+
if (A[i, 2] != 0) {
110+
y <- (b[i] - A[i, 1] * x) / A[i, 2]
111+
ylim <- range(c(ylim, y))
112+
}
113+
}
114+
}
115+
116+
labels <- if (isTRUE(labels)) {
117+
showEqn(A, b, vars, simplify=TRUE)
118+
}
119+
120+
for (i in 1:neq) {
121+
if (i == 1) plot(xlim, ylim, type="n", xlab = vars[1], ylab = vars[2], xlim = xlim, ylim = ylim, ...)
122+
123+
if (A[i, 2] == 0) {
124+
abline(v = b[i] / A[i, 1], col = col[i], lwd = lwd[i], lty = lty[i])
125+
y <- ylim
126+
}
127+
else {
128+
# calculate y values for current equation
129+
y <- (b[i] - A[i, 1] * x) / A[i, 2]
130+
lines(x, y, col = col[i], type = 'l', lwd = lwd[i], lty = lty[i])
131+
}
132+
133+
if (!is.null(labels)) {
134+
xl <- if(A[i, 2] == 0) b[i] else x[1]
135+
label <- parse(text=sub("=", "==", labels[i]))
136+
text(xl, y[1], label, col=col[i], pos=4)
137+
}
138+
}
139+
140+
if (axes) abline(h=0, v=0, col="gray")
141+
142+
if (!isFALSE(solution)) {
143+
if (is.list(solution)) {
144+
solution$cex <- solution$cex %||% 1.5
145+
solution$pch <- solution$cex %||% 16
146+
}
147+
points(intersections, cex = solution$cex, pch = solution$pch)
148+
}
149+
150+
}
151+
152+
153+
# plotEqn <- function(A, b, vars, xlim=c(-4, 4), ylim,
154+
# col=1:nrow(A), lwd=2, lty=1,
155+
# axes=TRUE, labels=TRUE,
156+
# solution=TRUE
157+
# ) {
158+
# if (!is.numeric(A) || !is.matrix(A)) stop("A must be a numeric matrix")
159+
# if (missing(b)) {
160+
# b <- A[,ncol(A)] # assume last column of Ab
161+
# A <- A[,-ncol(A)] # remove b from A
162+
# }
163+
# if (ncol(A) != 2) stop("plotEqn only handles two-variable equations. Use plotEqn3d for three-variable equations.")
164+
# if (missing(vars)) vars <- c(expression(x[1]), expression(x[2])) # paste0("x", 1:ncol(A))
165+
#
166+
# # set values for horizontal variable
167+
# x <- seq(xlim[1], xlim[2], length=10)
168+
#
169+
# neq <- nrow(A)
170+
# if (length(col) < neq) col <- rep_len(col, length.out=neq)
171+
# if (length(lwd) < neq) lwd <- rep_len(lwd, length.out=neq)
172+
# if (length(lty) < neq) lty <- rep_len(lty, length.out=neq)
173+
#
174+
# if (missing(ylim)) {
175+
# ylim <- xlim
176+
# for (i in 1:neq) {
177+
# if (A[i,2] != 0) {
178+
# y <- (b[i] - A[i,1] * x) / A[i,2]
179+
# ylim <- range(c(ylim, y))
180+
# }
181+
# }
182+
# }
183+
#
184+
# if (is.logical(labels) && labels) {
185+
# labels <- showEqn(A,b, vars, simplify=TRUE)
186+
# }
187+
# else labels=NULL
188+
#
189+
# for (i in 1:neq) {
190+
# if (i==1) plot(xlim, ylim, type="n", xlab = vars[1], ylab = vars[2], xlim = xlim, ylim = ylim)
191+
#
192+
# if (A[i,2] == 0) {
193+
# abline( v = b[i] / A[i,1], col = col[i], lwd = lwd[i], lty = lty[i] )
194+
# y <- ylim
195+
# }
196+
# else {
197+
# # calculate y values for current equation
198+
# y <- (b[i] - A[i,1] * x) / A[i,2]
199+
# lines( x, y, col = col[i], type = 'l', lwd = lwd[i], lty = lty[i] )
200+
# }
201+
#
202+
# if (!is.null(labels)) {
203+
# xl <- if(A[i,2] == 0) b[i] else x[1]
204+
# yl <- y[1]
205+
# label <- labels[i]
206+
# label <- parse(text=sub("=", "==", label))
207+
# text(xl, yl, label, col=col[i], pos=4)
208+
# }
209+
# }
210+
# if (axes) abline(h=0, v=0, col="gray")
211+
#
212+
# if (solution) {
213+
# for (i in 1:neq-1) {
214+
# for (j in i:neq) {
215+
# x <- try(solve(A[c(i,j),],b[c(i,j)]), silent=TRUE)
216+
# if (!inherits(x, "try-error")) points(x[1], x[2], cex=1.5)
217+
# }
218+
# }
219+
# }
220+
# }
221+
222+
223+
#' Plot Linear Equations in 3D
224+
#'
225+
#' Shows what matrices \eqn{A, b} look like as the system of linear equations, \eqn{A x = b} with three unknowns,
226+
#' x1, x2, and x3, by plotting a plane for each equation.
227+
228+
#' @param A either the matrix of coefficients of a system of linear equations, or the matrix \code{cbind(A,b)}
229+
#' The \code{A} matrix must have three columns.
230+
#' @param b if supplied, the vector of constants on the right hand side of the equations, of length matching
231+
#' the number of rows of \code{A}.
232+
#' @param vars a numeric or character vector of names of the variables.
233+
#' If supplied, the length must be equal to the number of unknowns in the equations.
234+
#' The default is \code{paste0("x", 1:ncol(A)}.
235+
#' @param xlim axis limits for the first variable
236+
#' @param ylim axis limits for the second variable
237+
#' @param zlim horizontal axis limits for the second variable; if missing, \code{zlim} is calculated from the
238+
#' range of the set of equations over the \code{xlim} and \code{ylim}
239+
#' @param col scalar or vector of colors for the lines, recycled as necessary
240+
#' @param alpha transparency applied to each plane
241+
#' @param labels logical, or a vector of character labels for the equations; not yet implemented.
242+
#' @param solution logical; should the solution point for all equations be marked (if possible)
243+
#' @param axes logical; whether to frame the plot with coordinate axes
244+
#' @param lit logical, specifying if lighting calculation should take place on geometry; see \code{\link[rgl]{rgl.material}}
245+
#'
246+
#' @return nothing; used for the side effect of making a plot
247+
#'
248+
#' @author Michael Friendly, John Fox
249+
#' @references Fox, J. and Friendly, M. (2016). "Visualizing Simultaneous Linear Equations, Geometric Vectors, and
250+
#' Least-Squares Regression with the matlib Package for R". \emph{useR Conference}, Stanford, CA, June 27 - June 30, 2016.
251+
#' @export
252+
#' @examples
253+
#' # three consistent equations in three unknowns
254+
#' A <- matrix(c(13, -4, 2, -4, 11, -2, 2, -2, 8), 3,3)
255+
#' b <- c(1,2,4)
256+
#' plotEqn3d(A,b)
257+
258+
plotEqn3d <- function( A, b, vars, xlim=c(-2,2), ylim=c(-2,2), zlim,
259+
col=2:(nrow(A)+1), alpha=0.9,
260+
labels=FALSE, solution=TRUE,
261+
axes=TRUE, lit=FALSE)
262+
{
263+
if (!is.numeric(A) || !is.matrix(A)) stop("A must be a numeric matrix")
264+
if (missing(b)) {
265+
b <- A[,ncol(A)] # assume last column of Ab
266+
A <- A[,-ncol(A)] # remove b from A
267+
}
268+
if (ncol(A) != 3) stop("plotEqn3d only handles three-variable equations")
269+
if (missing(vars)) vars <- paste0("x", 1:ncol(A))
270+
271+
neq <- nrow(A)
272+
# determine zlim if not specified
273+
if (missing(zlim)) {
274+
x <- xlim; y <- ylim
275+
zlim <- c(0, 0)
276+
for (i in 1:neq) {
277+
if (A[i,3] != 0) {
278+
z <- (b[i] - A[i,1] * x - A[i,2] * y) / A[i,3]
279+
zlim <- range(c(zlim, z))
280+
}
281+
}
282+
}
283+
284+
if (length(col) < neq) col <- rep_len(col, length.out=neq)
285+
286+
if (is.logical(labels) && labels) {
287+
# labels <- showEqn(A,b, vars)
288+
labels <- paste0("(", 1:neq, ")")
289+
}
290+
else labels=NULL
291+
292+
# rgl properties
293+
294+
depth_mask <- if (alpha < 1) TRUE else FALSE
295+
# Initialize the scene, no data plotted
296+
# Create some dummy data
297+
dat <- replicate(2, 1:3)
298+
rgl::plot3d(dat, type = 'n', xlim = xlim, ylim = ylim, zlim = c(-3, 3),
299+
xlab = vars[1], ylab = vars[2], zlab = vars[3],
300+
axes=axes)
301+
# Add planes
302+
rgl::planes3d(A[,1], A[,2], A[,3], -b,
303+
col=col, alpha=alpha, lit=lit, depth_mask=depth_mask)
304+
305+
# show the solution??
306+
if (solution) {
307+
x <- try(solve(A,b), silent=TRUE)
308+
if (!inherits(x, "try-error")) rgl::spheres3d(solve(A,b), radius=0.2)
309+
}
310+
311+
# if (!is.null(labels)) {
312+
# for (i in 1:neq) {
313+
# xl <- xlim[1]
314+
# yl <- ylim[1]
315+
# zl <- if (A[i,3] != 0) min((b[i] - A[i,1] * xl - A[i,2] * yl) / A[i,3]) else 0
316+
# rgl::text3d(xl, yl, zl, labels[i])
317+
# }
318+
# }
319+
}

‎dev/test-vectors3d.R

+3-2
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,8 @@
55
library(rgl)
66
library(matlib)
77

8-
# want to use labels that are expressions:
8+
# want to use labels that are expressions.
9+
# text3d() now allows this, using `plotmath3d()` for expressions
910
labs <- c(expression(x[1]), "y", expression(x[2]))
1011
is.expression(labs)
1112

@@ -16,7 +17,7 @@ E <- diag(3)
1617
rownames(E) <- c(expression(x[1]), "y", expression(x[2]))
1718
vectors3d(E, lwd=2)
1819
vectors3d(c(1, 1, 1),
19-
labels=c("",expression(hat(y))), color="red",
20+
labels=c(expression(hat(y))), color="red",
2021
lwd=3)
2122
vectors3d(c(1, 1, 0),
2223
labels=c("", "x+y"),

‎man/vectors3d.Rd

+2
Some generated files are not rendered by default. Learn more about customizing how changed files appear on GitHub.

0 commit comments

Comments
 (0)