Messing around with R

hello.r

# hello.r

# Windows one-time setup (requires choco):
#
#     choco install R
#
# Run this from a shell:
#
#     "/c/Program Files/R/R-4.1.2/bin/Rscript.exe" hello.r
#
# Or paste it into the online interpreter at:
#
#     https://rdrr.io/snippets/
#

#===============================================================================

# Seed the RNG (optional, only for repeatability)
set.seed(0)
#set.seed(1)

## \n is literal for print()!
#print("Hello world\n")

cat("\nHello world\n\n")

# Rectangular bounds of pointcloud distribution
dxh <- runif(1, 5.0 , 20.0)
dyh <- runif(1, 0.5 ,  2.0)
dzh <- runif(1, 0.05,  0.2)

# Rotation angles (radians)
t1 <- runif(1, 0, 2 * pi)
t2 <- runif(1, 0, 2 * pi)
t3 <- runif(1, 0, 2 * pi)

#===============================================================================
#
## This is me, learning R
#
##a <- array(1:6, c(3,2))
#
## Make a random 3x2 matrix.  Note total size vs nrows (no explicit cols)
#a <- matrix(runif(6), nrow=3)
#
#x <- array(1:2)
#
##dim(a)
##dim(x)
#
#cat("a =\n")
#a
#cat("\nx =\n")
#x
#
#b <- a %*% x
#
#cat("\na*x =\n")
#b
#===============================================================================

# Make np random points, uniformly distributed from -dxh to +dxh in x
# direction, -dyh to +dyh in y, ... .  Points are saved in matrix a.

# As np approaches \infty, eigenvalues/vectors of PCA should approach d[xyz]h
# and rotation matrix vectors (columns).  Standard deviations will depend on the
# distribution, e.g. normal or uniform

np <- 1000
nd <- 3

cat("dxh, dyh, dzh =\n")
dxh
dyh
dzh

## Setting a by using temp vecs
#x <- runif(np, -dxh, dxh)
#y <- runif(np, -dyh, dyh)
#z <- runif(np, -dzh, dzh)
#
## Compose a by catting columns
#a <- matrix(c(x, y, z), nrow = np)

# Skip the temp vecs

## Uniform distribution
#a <- matrix(c(runif(np, -dxh, dxh),
#              runif(np, -dyh, dyh),
#              runif(np, -dzh, dzh)), nrow = np)

# Normal distribution
a <- matrix(c(rnorm(np, -dxh, dxh),
              rnorm(np, -dyh, dyh),
              rnorm(np, -dzh, dzh)), nrow = np)

cat("dim(a) =\n")
dim(a)
#a

# Rotation matrices about each axis
r1 <- matrix(c(cos(t1), -sin(t1), 0,
               sin(t1),  cos(t1), 0,
                     0,        0, 1), nrow = nd)
#r1

r2 <- matrix(c(cos(t2), 0, -sin(t2),
                     0, 1,        0,
               sin(t2), 0,  cos(t2)), nrow = nd)
#r2

r3 <- matrix(c(1,       0,        0,
               0, cos(t3), -sin(t3),
               0, sin(t3),  cos(t3)), nrow = nd)
#r3

# Total rotation matrix.  Yes, that's really the operator they picked
r <- r1 %*% r2 %*% r3

# Transpose displays better w/ PCA output from prcomp().  Results may differ by
# a sign
t(r)

## This should be skew-symmetric for all rotation matrices r
#r - t(r)

# Rotate points
a <- a %*% r

#a
#t(a) %*% a

# Open RGui and run this script to see plot

#plot3d(a[,1], a[,2], a[,3])
plot(a[,1], a[,2])

# Run PCA on rotated data to reverse engineer the rotation of the original
# axis-aligned pointcloud, and the standard deviations along each dimension.
#
# This creates an object which can be assigned, but without assigning it just
# cats the results to the console.
#
prcomp(a)

#===============================================================================