stopifnot(require(FAiR)) # needs version 0.5-0 or greater
set.seed(12345)
n <- 10 # number of manifest variables
r <- 4  # number of inputs to data-generating process
if( r >= 0.5 * (2 * n + 1 - sqrt(8 *n + 1)) ) {
	print("Neither Theta nor functions thereof are identified")
}

# draw random communalities
h2 <- runif(n, min = .05, max = .95)

# draw random eigenvalues for Pi
pi <- rsimplex(r) * n / r
pi <- c(sort(pi, TRUE), rep(0, n - r))

# make random covariance matrix among manifest variables (or use your own data here)
man <- make_manifest(pi, communalities = h2)

# make restrictions objects
res.C <- make_restrictions(man, discrepancy = "THEIL", scaling = "c")
res.U <- make_restrictions(man, discrepancy = "THEIL", scaling = "u")

# find optimum
EFA.C <- Factanal(man, res.C, wait.generations = 20) # optimum wrt pi(Theta)
EFA.U <- Factanal(man, res.U, wait.generations = 20) # optimum wrt delta(Theta)

# check answer for eigenvalues
round( cbind(truth = pi, # true population eigenvalues
	C = eigen(cov2cor(fitted(EFA.C)), TRUE, TRUE)$values,
	U = eigen(cov2cor(fitted(EFA.U)), TRUE, TRUE)$values), 3)

# check answer for Theta
round( cbind(truth = 1 - h2, # true error variances
	C = uniquenesses(EFA.C),
	U = uniquenesses(EFA.U)), 3 )

# steak sauce