info   0
) Macros for doing ULS and GLS factor analysis estimation
)
) mlcrit()  Macro to compute ML criterion function and gradient
) glscrit() Macro to compute GLS criterion function and gradient
) ulscrit() Macro to compute ULS criterion function and gradient
) facanal() Macro to use optimization macro dfp(), bfs() or broyden()
)    to find ULS, GLS or ML estimates of uniquenesses and loadings.
) mlfactor() Macro to optimization macro to find ML estimates of 
)    uniquenesses and loadings
) ulsfactor(s or r, m [[, quiet:T,silent:T,start:psi0,uselogs:T],\
)       maxit:nmax, minit:nmin, print:T,crit:crvec])
)   does ULS factor extraction by non-linear regression
) ulsresids(psi,0,vector(s or r),m)
)   Macro used by ulsfactor
) glsfactor(s or r, m [[, quiet:T,silent:T,start:psi0,uselogs:T],\
)       maxit:nmax, minit:nmin, print:T,crit:crvec])
)   does GLS factor extraction by non-linear regression
) glsresids(psi,vector(s or r),vector(s or r),m)
)   Macro used by glsfactor
)
) The following macros essentially duplicate the factor analysis macros
) in mulvar.mac which were updated versions of macros in MacAnova.mac
)
)   goodfit(s,lhat,psihat)
)   stepuls(s,psi,m [,print:T])
)   stepml(s,psi,m [,print:T])
)   stepgls(s,psi,m [,print:T])
)
))File version 011120

mlcrit       MACRO DOLLARS
) Macro to compute ML criterion function and its gradient (computed from
) differences)
) The parameters x are log(psi_sub_i), i = 1,...,p
) This is to be used with optimization macros dfp(), bfs() and broyden()
) Usage
)  result <- mlcrit(logpsi, k, structure(s, nfact [, del]))
)   logpsi       REAL vector of log uniquenesses
)   k            integer scalar to control what is calculated
)                k = 0, return criterion (scalar) (MISSING when s not P.D.)
)                k = 1, return gradient (vector) (rep(MISSING,nrows(s))
)                  when s not P.D.)
)                k = -2, return structure(loadings, eigenvals)
)   s            REAL symmetric covariance or correlation matrix
)   nfact > 0    integer number of factors
)   del >= 0     Small REAL scalar [1e-5].  Gradient is computed as 
)     gradient[i] = (critfn(logpsi + del*(run(p)==i)) - critfn(logpsi))/del
)                del = 0 results in analytic gradient being used
)   When s not positive definite, MISSING (k = 0) or
)   rep(MISSING,nrows(s)) (k = 1) is returned
) No argument checking is done
)) 011120 First draft
)) 011124 del = 0 signals use of analytic gradient
)) Version 011124, C. Bingham, kb@stat.umn.edu
@x <- $1
@n <- length(@x)
_PSIML <- exp(@x)

@param <- $3
@s <- @param[1]
@J <- run(@param[2])
@result <- if (($2) < 0) {
	if (($2) == -2) {
	  # return structure(loadings, eigenvals)
		@stuff <- releigen(@s, dmat(_PSIML))
		@gamma <- @stuff$values[@J]-1
		@crit <- NULL
		@loadings <- if(min(@gamma) < 0){
			NULL
		} else {
			_PSIML * @stuff$vectors[,@J] * sqrt(@gamma')
		}
		structure(loadings:delete(@loadings,return:T),\
				  eigenvals:delete(@stuff,return:T)$values)
	} else {
		NULL
	}
} else {
	@stuff <- releigen(@s, dmat(_PSIML))
	@gamma <- @stuff$values[-@J]
	if (min(@gamma) <= 0) {
		delete(@stuff)
		if (($2) == 0) {
			?
		} else {
			rep(?, @n)
		}
	} else {
		@crit <- sum(@gamma - 1 - log(@gamma))
		if (($2) == 0){
			@crit
		} else {
			@del <- if (ncomps(@param) > 2) {
				@param[3]
			} else {
				1e-5
			}
			if (@del > 0) {
				@gradient <- @x
				for (@i,1,@n){
					@tmp <- @x
					@tmp[@i] <- @tmp[@i] + @del
					@gamma <- releigenvals(@s, dmat(exp(@tmp)))[-@J]
					@gradient[@i] <- if (min(@gamma) <= 0) {
						?
					} else {
						(sum(@gamma - 1 - log(@gamma)) - @crit)/@del
					}
				}
				delete(@tmp,@i)
			} else {
				@omega <- sqrt(_PSIML)*@stuff$vectors[,-@J]
				@gradient <- vector(sum((1 - @gamma)*(@omega'^2)))
			}
			
			delete(@del,@stuff,@crit)
			delete(@gradient,return:T)
		}
	}
}

delete(@x,@s,@gamma,@crit,@J,@param,@n,silent:T)
delete(@result,return:T)
%mlcrit%

glscrit       MACRO DOLLARS
) Macro to compute GLS criterion function and its gradient (computed from
) differences)
) The parameters x are log(psi_sub_i), i = 1,...,p
) This is to be used with optimization macros dfp(), bfs() and broyden()
) Usage
)  result <- glscrit(logpsi, k, structure(s, nfact [, del]))
)   logpsi       REAL vector of log uniquenesses
)   k            integer scalar to control what is calculated
)                k = 0, return criterion (scalar) (MISSING when s not P.D.)
)                k = 1, return gradient (vector) (rep(MISSING,nrows(s))
)                  when s not P.D.)
)                k = -2, return structure(loadings, eigenvals)
)   s            REAL symmetric covariance or correlation matrix
)   nfact > 0    integer number of factors
)   del >= 0     Small REAL scalar [1e-5].  Gradient is computed as 
)     gradient[i] = (critfn(logpsi + del*(run(p)==i)) - critfn(logpsi))/del
)                del = 0 results in analytic gradient being used
)   When s not positive definite, MISSING (k = 0) or
)   rep(MISSING,nrows(s)) (k = 1) is returned
) No argument checking is done
)) Version 011120, C. Bingham, kb@stat.umn.edu
)) 011124 del = 0 signals use of analytic gradient
@x <- $1
@n <- length(@x)
_PSIGLS <- exp(@x)

@param <- $3
@s <- @param[1]
@J <- run(@param[2])
@result <- if (($2) < 0) {
	if (($2) == -2) {
		# return structure(loadings, eigenvalues)
		@stuff <- releigen(@s, dmat(_PSIGLS))
		@gamma <- @stuff$values[@J]-1
		@crit <- NULL
		@loadings <- if(min(@gamma) < 0){
			NULL
		} else {
			_PSIGLS * @stuff$vectors[,@J] * sqrt(@gamma')
		}
		structure(loadings:delete(@loadings,return:T),\
				  eigenvals:delete(@stuff,return:T)$values)
	} else {
		NULL
	}
} else {
	@stuff <- releigen(@s, dmat(_PSIGLS))
	@gamma <- @stuff$values[-@J]
	if (min(@gamma) <= 0) {
		delete(@stuff)
		if (($2) == 0) {
			?
		} else {
			rep(?, @n)
		}
	} else {
		@crit <- .5*sum((1 - 1/@gamma)^2)
		@result <- if (($2) == 0){
			@crit
		} else {
			@del <- if (ncomps(@param) > 2) {
				@param[3]
			} else {
				1e-5
			}
			if (@del > 0) {
				@gradient <- @x
				for (@i,1,@n){
					@tmp <- @x
					@tmp[@i] <- @tmp[@i] + @del
					@gamma <- releigenvals(@s, dmat(exp(@tmp)))[-@J]
					@gradient[@i] <- if (min(@gamma) <= 0) {
						?
					} else {
						@gradient[@i] <- (.5*sum((1 - 1/@gamma)^2) - @crit)/@del
					}
				}
				delete(@tmp,@i)
			} else {
				@omega <- sqrt(_PSIGLS)*@stuff$vectors[,-@J]
				@gradient <- vector(sum((1 - @gamma)*(@omega'^2)/@gamma^2))
			}
			delete(@del,@stuff)
			delete(@gradient,return:T)
		}
	}
}
delete(@x,@s,@gamma,@crit,@J,@n,silent:T)
delete(@result,return:T)
%glscrit%

ulscrit       MACRO DOLLARS
) Macro to compute ULS criterion function and its gradient (computed from
) differences)
) The parameters x are log(psi_sub_i), i = 1,...,p
) This is to be used with optimization macros dfp(), bfs() and broyden()
) Usage
)  result <- ulscrit(logpsi, k, structure(s, nfact [, del]))
)   logpsi       REAL vector of log uniquenesses
)   k            integer scalar to control what is calculated
)                k = 0, return criterion (scalar) (MISSING when s not P.D.)
)                k = 1, return gradient (vector) (rep(MISSING,nrows(s))
)                  when s not P.D.)
)                k = -2, return structure(loadings, eigenvals)
)   s            REAL symmetric covariance or correlation matrix
)   nfact > 0    integer number of factors
)   del >= 0     Small REAL scalar [1e-5].  Gradient is computed as 
)     gradient[i] = (critfn(logpsi + del*(run(p)==i)) - critfn(logpsi))/del
)                Currently del = 0 gets the default.  At some point del = 0
)                will signal the use of analytic derivatives
)   When s not positive definite, MISSING (k = 0) or
)   rep(MISSING,nrows(s)) (k = 1) is returned
) No argument checking is done
)) Version 011120, C. Bingham, kb@stat.umn.edu
@x <- $1
@n <- length(@x)
_PSIULS <- exp(@x)

@param <- $3
@s <- @param[1]
@J <- run(@param[2])
@result <- if (($2) < 0) {
	if (($2) == -2) {
		# return structure(loadings, eigenvals)
		@crit <- NULL
		@stuff <- eigen(@s - dmat(_PSIULS))
		@delta <- @stuff$values[@J]
		@loadings <- if(min(@delta) < 0){
			NULL
		} else {
			@stuff$vectors[,@J] * sqrt(@delta')
		}
		structure(loadings:delete(@loadings,return:T),\
				  eigenvals:delete(@stuff,return:T)$values)
	} else {
		NULL
	}
} else {
	@stuff <- eigen(@s - dmat(_PSIULS))
	@delta <- @stuff$values[-@J]
	@crit <- .5*sum(@delta^2)
	if (($2) == 0){
		@crit
	} else {
		@del <- if (ncomps(@param) > 2) {
			@param[3]
		} else {
			1e-5
		}
		if (@del > 0) {
			delete(@stuff)
			@gradient <- @x
			for (@i,1,@n){
				@tmp <- @x
				@tmp[@i] <- @tmp[@i] + @del
				@delta <- eigenvals(@s - dmat(exp(@tmp)))[-@J]
				@gradient[@i] <- (.5*sum(@delta^2) - @crit)/@del
			}
			delete(@tmp,@i,@del)
		} else {
			@omega <- delete(@stuff,return:T)$vectors[,-@J]
			@gradient <-\
			   -_PSIULS*vector(sum(@delta*(delete(@omega,return:T)'^2)))
		}
		
		delete(@gradient,return:T)
	}
}

delete(@x,@s,@delta,@crit,@J,@n,silent:T)
delete(@result,return:T)
%ulscrit%

facanal   MACRO   DOLLARS
) Macro to use optimization macro dfp(), bfs() or broyden()
) to find ULS, GLS or ML estimates of uniquenesses and loadings.
) It uses macro dfp(), bfs() or broyden() to minimize the criterion as
) a function of log(psi), psi = uniquenesses.
)
) It requires macros ulscrit(), glscrit() or mlcrit() to compute
) the objective function, gradient vector and loading matrix.
) 
) Usage:
)  facanal(s, m [, method:Meth] [,rotate:Rotmeth], [,start:psi0] \
)           [,maxit:nmax] [,minit:nmin] [,crit:vector(nsig1,nsig2,dg)]\
)           [,minimizer:M] [,gold:ngold] [,quiet:T] [,silent:T] \
)           [,recordwhen:d1] [,printwhen:d2])
)   s          REAL p x p positive definite covariance or correlation
)              matrix
)   m          0 < integer < (2*p + 1 - sqrt(8*@p+1))/2
)   Meth       one of "uls", "gls", "ml" or "mle" ["mle"]
)   Rotmeth    one of "varimax", "quartimax", "equimax" ["none"]
)   psi0       Positive length p vector of starting values for psi
)   nmax       integer > 0, maximum number of iterations permitted [30]
)   nmin       integer >= 0, minimum number of iterations [0]
)   nsig1      target number of significant digits in log uniquenesses [5]
)   nsig2      target number of significant digits in criterion [8]
)   dg         target upper limit for ||gradient|| [-1]
)   M          one of "dfp", "bfs", or "broyden" ["bfs"]
)   ngold      integer >= 0, # of steps in golden mean linear search [1]
)   quiet:T    suppresses printing of results
)   silent:T   suppresses printing of results and warning messages
)   d1         integer >= 0. When d1 > 0, partial results are printed at
)              iterations d1, 2*d1, ...  [0]; not affected by silent:T
)   d2         integer >= 0. When d2 > 0, partial results are recorded at
)              iterations d2, 2*d2, ...  [0]
)
)   nsig1, nsig2 and dg specify convergence critera; any <= 0 are ignored
)   ngold is ignored if M is "broyden"
) The result is an "invisible" structure
) structure(psihat:psi,loadings:l,criterion:crit,eigenvals:vals,\
)           gradient:g, method:Meth, rotation:Rotmeth,iter:n, status:k)
)   
)  psi         REAL vector of estimated uniquenesses
)  l           REAL p by m matrix of loadings
)  crit        Minimized criterion
)  vals        eigenvalues s relative to psi
)  g           REAL gradient vector at optimum
)  n > 0       integer = number of iterations
)  k >= 0      integer specifying which of the convergence criterion
)              signalled convergence; k == 0 means not converged
)
) When d2 > 0, side effect structure DFPRECORD, BNFRECORD or
) BROYDNRECORD with components funvals, xvals and gradients is created.
) On each call, ulscrit(), glscrit() or mlcrit() saves a copy of psi
) in invisible variable _PSIULS, _PSIGLS or _PSIML.
) 
)) Version of 011117, written by C. Bingham, kb@stat.umn.edu
)) 011125 added keyword 'golden' with default 1
))        added keyword 'rotation' with default "none"
))        changed default for delta to 0 so as to use analytic gradient
))        modified printed output
))        result now includes method and rotation
# $S(s or r, m [,start:psi0,maxit:m,silent:T,quiet:T])
@s <- argvalue($1,"variance or correlation","square real nonmissing")
@m <- argvalue($2,"number of factors","positive count")
@keys <- if ($k > 0) {
	structure($K)
} else {
	structure(NoTaKey:NULL)
}
@method <- keyvalue(@keys,"method","string",default:"ml")
@maxit <- keyvalue(@keys, "maxit*","positive count",default:30)
@minit <- keyvalue(@keys, "minit*","count",default:0)
@delta <- keyvalue(@keys,"delta","nonnegative number",default:0)
@silent <- keyvalue(@keys, "silent", "TF", default:F)
@quiet <- keyvalue(@keys, "quiet", "TF", default:@silent)
@record <- keyvalue(@keys, "recordwhen", "count", default:0)
@print <- keyvalue(@keys, "printwhen", "count", default:0)
@minimizer <- keyvalue(@keys,"minimizer","string",default:"bfs")
@crit <- keyvalue(@keys,"crit","real nonmissing vector",default:vector(5,8,-1))
@gold <- keyvalue(@keys,"gold*","count",default:1)
@rotmethod <- keyvalue(@keys,"rot*","string",default:"none")
@kaiser <- keyvalue(@keys,"kais*","TF",default:T)
delete(@keys)

@imeth <- min(match(@method,vector("uls","gls","ml","mle"),0),3)
if (@imeth == 0) {
	error(paste("unrecognized method \"",@method,"\"",sep:""))
}
if (sum(match(vector("var*","quart*","equi*","none"),\
			  @rotmethod,0,exact:F)) != 1) {
	error(paste(@rotmethod,\
				"is not one of \"varimax\",\"quartimax\" or \"equimax\""))
}

if (@minimizer == "broyden") {
	@gold <- 0
}

@p <- nrows(@s)
@kind <- if(sum(diag(@s) == 1) == @p){
	"correlation"
}else{
	"covariance"
}
if (min(eigenvals(@s)) <= 0){
	error(paste(@kind,"matrix is not positive definite"))
}
if (F && !ismacro(<<@minimizer>>)){
	if (@minimizer == "dfp") {
		getmacros(dfp,quiet:T,printname:F)
	} elseif (@minimizer == "bfs") {
		getmacros(bfs,quiet:T,printname:F)
	} elseif (@minimizer == "broyden") { 
		getmacros(broyden,quiet:T,printname:F)
	} else {
		error(paste("unrecognized 'minimizer' \"",@minimizer,"\"",sep:""))
	}
} elseif (!ismacro("minimizer")){
		getmacros(minimizer,quiet:T,printname:F)
}

@critfn <- vector("ulscrit","glscrit","mlcrit")[@imeth]
if (!ismacro(<<@critfn>>)) {
	if (@imeth == 1){
		getmacros(ulscrit,quiet:T,printname:F)
	} elseif (@imeth == 2){
		getmacros(glscrit,quiet:T,printname:F)
	} else {
		getmacros(mlcrit,quiet:T,printname:F)
	}
}

@x0 <- log(keyvalue($K,"start","positive vector",\
		   default:1/diag(solve(@s)))) # starting value

if (length(@x0) != @p){
	error("Wrong number of starting uniquenesses")
}

@param <- structure(@s, @m, @delta)
if (F) {
	@result <- <<@minimizer>>(@x0,<<@critfn>>,@param,\
		maxit:@maxit, minit:@minit, crit:@crit,gold:@gold,\
		recordwhen:@record, printwhen:@print)
} else {
	@result <- minimizer(@x0,<<@critfn>>,@param,method:@minimizer,\
		maxit:@maxit, minit:@minit, crit:@crit,gold:@gold,\
		recordwhen:@record, printwhen:@print,silent:@silent)
}

 # The result is structure(x:xmin, f:minVal, gradient:gradient,\
 # h:invhessian, iterations:niter)

@iter <- @result$iterations
@status <- @result$status
if (@status < 0) {
	error(paste("criterion or gradient not defined on iteration", @iter))
}

@x <- @result$x
PSI <- vector(exp(@x))
CRITERION <- @result$f
if(@method == "ml" || @method == "mle") {
	@tmp <- releigen(@s,dmat(PSI))$values[-run(@m)]
	CRITERION <- sum(@tmp-1-log(@tmp))
} 
GRADIENT <- vector(@result$gradient)

@result <- if (@imeth == 1) {
	ulscrit(@x, -2, @param)
} elseif (@imeth == 2){
	glscrit(@x, -2, @param)
} else {
	mlcrit(@x, -2, @param)
}
LOADINGS <- if (@rotmethod == "none") {
	@result$loadings
} else {
	@H <- if (!@kaiser) {
		1
	} else {
		vector(sqrt(sum(@result$loadings'^2)))
	}
	@H*rotation(@result$loadings/@H,method:@rotmethod,reorder:T)
}

EIGENVALUES <- delete(@result,return:T)$eigenvals
delete(@x,@param,@kind)

if (haslabels(@s)){
	setlabels(PSI, getlabels(@s,1))
	if (!isnull(LOADINGS)) {
		setlabels(LOADINGS, structure(getlabels(@s,1),"Factor "))
	}
}

if (!@silent) {
	if (2*@m >= 2*@p + 1 - sqrt(8*@p+1)){
		print(paste("WARNING: With m = ",@m," and p = ",@p,\
					", m > (2*p + 1 - sqrt(8*p+1))/2",sep:""))
	}
	if (isnull(LOADINGS)){
		print(paste("WARNING: eigenvalue",@m,"of Psi^(-1) %*% (R or S) < 1"))
	}
	if (@status <= 0){
		print(paste("WARNING: no convergence in",@iter,"iterations"))
	}
	if (!@quiet) {
		if(@status > 0){
			print(paste("Convergence in",@iter,"iterations by criterion",@status))
		}
		@loadname <- if (@rotmethod != "none") {
			paste(@rotmethod,"rotated estimated loadings")
		} else {
			"unrotated estimated loadings"
		}

		print(PSI,name:"estimated uniquenesses")
		print(LOADINGS, name:delete(@loadname,return:T))
		print(CRITERION, name:paste("minimized",@method,"criterion"))
	}
}

delete(@quiet,@silent,@p,@m,@delta)
@result <- structure(psihat:PSI,loadings:LOADINGS,criterion:CRITERION,\
	eigenvals:EIGENVALUES, gradient:GRADIENT,\
	method:delete(@method,return:T),rotation:delete(@rotmethod,return:T),\
	iter:delete(@iter,return:T),status:delete(@status,return:T))
delete(@result,return:T,invisible:T)
%facanal%

====> mlfactor <====
mlfactor   MACRO   DOLLARS
) Macro to use optimization macro dfp(), bfs() or broyden()
) to find ML estimates of uniquenesses and loadings.
)
) It requires macro mlcrit() to compute the objective function and
) gradient vector.
) Usage:
)  mlfactor(s, m [, start:psi0, maxit:nmax, minit:nmin])
)   s          REAL p x p positive definite covariance or correlation
)              matrix
)   m          0 < integer < (2*p + 1 - sqrt(8*@p+1))/2
)   psi0       Positive length p vector of starting values for psi
)   nmax       integer > 0, maximum number of iterations permitted [30]
)
) The result is an "invisible" structure
) structure(psihat:psi,loadings:l,criterion:crit,eigenvals:vals,\
)           gradient:g,iter:n)
)   
)  psi       REAL vector of estimated uniquenesses
)  l         REAL p by m matrix of loadings
)  crit      Minimized criterion
)  g         REAL gradient vector at optimum
)  vals      eigenvalues s relative to psi
)
) On each call, mlfactor() saves a copy of psi in invisible
) variable _PSIML.  Type 'print(_PSIML) to see this value.
)) Version of 011112, written by C. Bingham, kb@stat.umn.edu
)) It uses macro dfp() to minimize the ML criterion.
# $S(s or r, m [,start:psi0,maxit:m,silent:T,quiet:T])
@s <- argvalue($1,"variance or correlation","square real nonmissing")
@m <- argvalue($2,"number of factors","positive count")

@maxit <- keyvalue($K, "maxit*","positive integer",default:30)
@eps <- keyvalue($K, "eps*","positive number",default:1e-7)
@delta <- keyvalue($K,"delta","positive number",default:1e-5)
@silent <- keyvalue($K, "silent", "TF", default:F)
@quiet <- keyvalue($K, "quiet", "TF", default:T)
@record <- keyvalue($K, "recordwhen", "count", default:0)
@print <- keyvalue($K, "printwhen", "count", default:0)
@minimizer <- keyvalue($K,"minimizer","string",default:"dfp")
@p <- nrows(@s)
@kind <- if(sum(diag(@s) == 1) == @p){
	"correlation"
}else{
	"covariance"
}
if (min(eigenvals(@s)) <= 0){
	error(paste(@kind,"matrix is not positive definite"))
}
if (2*@m >= 2*@p + 1 - sqrt(8*@p+1)){
	error(paste("With m = ",@m," and p = ",@p,\
		", m > (2*p + 1 - sqrt(8*p+1))/2",sep:""))
}

if (!ismacro(<<@minimizer>>)){
	if (@minimizer == "dfp") {
		getmacros(dfp,quiet:T,printname:F)
	} elseif (@minimizer == "bfs") {
		getmacros(bfs,quiet:T,printname:F)
	} elseif (@minimizer == "broyden") { 
		getmacros(broyden,quiet:T,printname:F)
	} else {
		error(paste("unrecognized value",@minimizer,"for 'minimizer'"))
	}
}

if (!ismacro(mlcrit)){
	getmacros(mlcrit,quiet:T,printname:F)
}
@x0 <- log(keyvalue($K,"start","positive vector",\
	default:1/diag(solve(@s)))) # starting value

if (length(@x0) != @p){
	error("Wrong number of starting uniqueness")
}

@result <- <<@minimizer>>(@x0,mlcrit,structure(@s,@m,@delta),\
	maxit:@maxit, eps:@eps, recordwhen:@record, printwhen:@print)
 # The result is structure(x:xmin, f:minVal, gradient:gradient,\
 # h:invhessian, iterations:niter)

PSI <- exp(@result$x)

if (haslabels(@s)){
	setlabels(PSI, getlabels(@s,1))
}

@iter <- @result$iterations
CRITERION <- @crit <- @result$f
@gradient <- @result$gradient
delete(@result)

@eigs <- releigen(@s, dmat(PSI))
@J <- run(@m)
@vals <- @eigs$values[@J]
LOADINGS <- if(min(@vals) < 1){
	NULL
} else {
	PSI * @eigs$vectors[,@J] * sqrt(@vals-1)'
}

if (!@silent){
	if (isnull(LOADINGS)){
		print(paste("WARNING: eigenvalue",@m,"of Psi^(-1) %*% (R or S) < 1"))

	}
	if (@iter > @maxit){
		print(paste("WARNING: no convergence in",@iter,"iterations"))
	}
}
if (!@quiet && !@silent){
	if(@iter <= @maxit){
		print(paste("Convergence in",@iter,"iterations"))
	}	
	print(psihat:PSI,loadings:LOADINGS,criterion:CRITERION,eigenvals:@eigs$values)
}

delete(@quiet,@silent,@p,@m,@J,@delta)
@result <- structure(psihat:PSI,loadings:LOADINGS,criterion:CRITERION,\
	eigenvals:@eigs$values,gradient:delete(@gradient, return:T),\
	iter:delete(@iter,return:T))
delete(@eigs)
delete(@result,return:T,invisible:T)
%mlfactor%

====> ulsfactor <====
ulsfactor   MACRO   DOLLARS
) Macro to use nonlinear least squares macro levmar to find
) ULS estimates of uniquenesses and loadings.  It uses macro ulsresids
) to compute residuals used by levmar.
)
) Usage:
)	ulsfactor(s, m [, quiet:T,silent:T,start:psi0,uselogs:T],\
)       maxit:nmax, minit:nmin, print:T,crit:crvec])
)   s              positive definite covariance or correlation matrix
)   m              positive integer < (2*p + 1 - sqrt(8*@p+1))/2
)   nmax           nonnegative integer, maximum number of iterations
)                  permitted [30]
)   nmin           nonnegative integer, minimum of iterations required [0]
)   crvec          vector(numsig, nsigsq, delta), 3 criteria for
)                  convergence (a negative criterion is not used)
)   psi0           positive REAL vector of starting uniquenesses
)   With print:T, something is printed on each interation
)
) There is no guarantee that min(psihat) > = 0.  In that case you can try
)   using uselogs:T to parametrize in terms of log(psi).  This ensures
)    min(psi) > 0 if it converges.
)   Howerever, if the minimum is actually outside the permissible
)    range, it will probably lead to a 'argument to solve() singular'
)    error message.
)
) The result is an "invisible" structure
) structure(psihat:psi,loadings:l,criterion:crit,\
)   eigenvals:vals,status:iconv,iter:n)
) where
)  psi               vector of estimated uniquenesses
)  l                 p by m matrix of loadings satisfying
)  crit              sum of squared differences of s from sigmathat
)  vals              vector of eigenvalues eigenvalues of s - psi
)  iconv = 0 or 4 not converged
)        = 1, converged with relative change in coeffs < relcon = 10^-numsig
)        = 2, converged with relative change in rss <= relssq = 10^-nsigsq
)        = 3, converged with ||gradient|| < delta
)        = 4, interation stopped; could not reduce rss.
)
) Macro ulsresids saves a copy of psi on every call in invisible variable
) _PSIULS.  Type 'print(_PSIULS) to see this value.
)
)) Version of 991114, written by C. Bingham, kb@stat.umn.edu
)) It uses macro levmar to minimize the sum of squared residuals.
# $S(s or r, m [,start:psi0,maxit:m,print:T,quiet:T])
@s <- argvalue($1,"variance or correlation","square real nonmissing")

@quiet <- keyvalue($K, "quiet","TF",default:F)
@silent <- keyvalue($K, "silent","TF",default:F)
@maxit <- keyvalue($K, "maxit*","positive integer",default:30)

__USELOGSULS <- keyvalue($K, "uselogs","TF",default:F)

@m <- argvalue($2,"number of factors","positive integer")
@p <- nrows(@s)
@kind <- if(sum(diag(@s) == 1) == @p){
	"correlation"
}else{
	"covariance"
}
if (min(eigenvals(@s)) <= 0){
	error(paste(@kind,"matrix is not positive definite"))
}
if (2*@m >= 2*@p + 1 - sqrt(8*@p+1)){
	error(paste("With m = ",@m," and p = ",@p,\
		", m > (2*p + 1 - sqrt(8*p+1))/2",sep:""))
}

if (!ismacro(levmar)){
	getmacros(levmar,quiet:T,printname:F)
}
if (!ismacro(ulsresids)){
	getmacros(ulsresids,quiet:T,printname:F)
}

@psi0 <- keyvalue($K,"start","positive vector",\
	default:1/diag(solve(@s))) # starting value
if (length(@psi0) != @p){
	error("Wrong number of starting uniqueness")
}
if (__USELOGSULS){
	@psi0 <- log(@psi0)
}
@result <- levmar(@psi0,0,vector(@s),@m,resid:ulsresids,$K)

@psihat <- if (__USELOGSULS){
	exp(@result$coefs)
}else{
	@result$coefs
}
@iconv <- @result$iconv
@iter <- @result$iter

if (@iconv == 0){
	print(paste("WARNING: failed to converge in",@maxit,"iterations"))
}
@eigs <- eigen(@s - dmat(@psihat))
@vals <- @eigs$values
@J <- run(@m)
@lhat <- if (@vals[@m] < 0){
	NULL
}else{
	@eigs$vectors[,@J] * sqrt(@vals[@J]')
}
delete(@eigs,@J)

if (!@silent){
	if (isnull(@lhat)){
		print(paste("WARNING: eigenvalue",@m,"of (R or S) - Psi is negative"))

	}
	if (min(@psihat) < 0){
		print("WARNING: psi not acceptable since min(psihat) < 0")
	}
	if (@iconv == 4){
		print("WARNING: could not reduce RSS on iteration",@iter)
	}elseif (@iconv <= 0){
		print(paste("WARNING: no convergence in",@iter,"iterations"))
	}
}
if (!@quiet && !@silent){
	if(@iconv > 0 && @iconv < 4){
		print(paste("Convergence in",@iter,"iterations by criterion",@iconv))
	}	
	print(psihat:@psihat,loadings:@lhat,criterion:@result$rss,eigenvals:@vals)
}

delete(@quiet,@silent,@p,@m,@psi0,__USELOGSULS)
@result <- structure(psihat:delete(@psihat,return:T),\
	loadings:delete(@lhat,return:T),criterion:delete(@result$rss,return:T),\
	eigenvals:delete(@vals,return:T),status:delete(@iconv,return:T),\
	iter:delete(@iter,return:T))
delete(@result,return:T,invisible:T)
%ulsfactor%

====> ulsresids <====
ulsresids    MACRO    DOLLARS OUTLINE
) Macro to be used with levmar to do ULS factor extraction
)
) levmar usage:
)   levmar(if(__USELOGSULS){log(psi)}else{psi},0,\
)     vector(r),m,resid:ulsresids)
)) 0000828
# $S(psi,0,vector(r),m)
_PSIULS <- if(__USELOGSULS){
	exp($1)
}else{
	$1
}
@psi <- dmat(_PSIULS)
@r <- $3
@p <- sqrt(length(@r))
@J <- run($4) #run(m)
@eigs <- eigen(matrix(@r,@p)  - @psi)
@l <- @eigs$vectors[,@J]
@res <- delete(@r,return:T) - vector((@l * @eigs$values[@J]') %C% @l + @psi)
delete(@psi,@p,@J,@eigs,@l)
delete(@res, return:T)
%ulsresids%

====> glsfactor <====
glsfactor   MACRO   DOLLARS
) Macro to use nonlinear least squares macro levmar to find
) GLS estimates of uniquenesses and loadings.  It uses macro glsresids
) to compute residuals used by levmar.
)
) Usage:
)	glsfactor(s, m [, quiet:T,silent:T,start:psi0,uselogs:T],\
)       maxit:nmax, minit:nmin, print:T,crit:crvec])
)   s: a positive definite covariance or correlation matrix
)   m: a positive integer < (2*p + 1 - sqrt(8*p+1))/2
)   nmax: nonnegative integer, maximum number of iterations permitted [30]
)   nmin: nonnegative integer, minimum of iterations required [0]
)   crvec: vector(numsig, nsigsq, delta), 3 criteria for convergence,
)    a negative criterion is not used
)   With print:T, something is printed on each interation
)
) There is no guarantee that min(psihat) > = 0.  In that case you can try
)   using uselogs:T to parametrize in terms of log(psi).  This ensures
)    min(psi) > 0 if it converges.
)   Howerever, but if the minimum is actually outside the permissible
)    range, it will probably lead to a 'argument to solve() singular'
)    error message.
)
) The result is an "invisible" structure
)  structure(psihat:psi,loadings:l,criterion:crit,\
)      eigenvals:vals,status:iconv,iter:n)
) where
)  psi is the vector of estimated uniquenesses
)  l is the p by m matrix of loadings satisfying
)  crit is sum of squared differences of s from sigmathat
)  vals are eigenvalues of s - psi
)  iconv = 0 or 4 not converged
)        = 1, converged with relative change in coeffs < relcon = 10^-numsig
)        = 2, converged with relative change in rss <= relssq = 10^-nsigsq
)        = 3, converged with ||gradient|| < delta
)        = 4, interation stopped; could not reduce rss.
)
) Macro glsresids saves a copy of psi on every call in invisible variable
) _PSIGLS.  Type 'print(_PSIGLS) to see this value.
)) Version of 991115, written by C. Bingham, kb@stat.umn.edu
)) It uses macro levmar to minimize the sum of squared residuals.
# $S(s or r, m [,start:psi0,maxit:m,print:T,quiet:T])
@s <- argvalue($1,"variance or correlation","square real nonmissing")

@quiet <- keyvalue($K, "quiet","TF",default:F)
@silent <- keyvalue($K, "silent","TF",default:F)
@maxit <- keyvalue($K, "maxit*","positive integer",default:30)

__USELOGSGLS <- keyvalue($K, "uselog*","TF",default:F)

@m <- argvalue($2,"number of factors","positive integer")
@p <- nrows(@s)
@kind <- if(sum(diag(@s) == 1) == @p){
	@what <- "r"
	"correlation"
}else{
	@what <- "s"
	"covariance"
}
@kind <- if(sum(diag(@s) == 1) == @p){
	"correlation"
}else{
	"covariance"
}
if (min(eigenvals(@s)) <= 0){
	error(paste(@kind,"matrix is not positive definite"))
}
if (2*@m >= 2*@p + 1 - sqrt(8*@p+1)){
	error(paste("With m = ",@m," and p = ",@p,\
		", m > (2*p + 1 - sqrt(8*p+1))/2",sep:""))
}

if (!ismacro(levmar)){
	getmacros(levmar,quiet:T,printname:F)
}
if (!ismacro(glsresids)){
	getmacros(glsresids,quiet:T,printname:F)
}

@psi0 <- keyvalue($K,"start","positive vector",\
	default:1/diag(solve(@s))) # starting value
if (length(@psi0) != @p){
	error("Wrong number of starting uniqueness")
}
if (__USELOGSGLS){
	@psi0 <- log(@psi0)
}
@result <- levmar(@psi0,vector(@s),vector(@s),@m,resid:glsresids,$K)

@psihat <- if (__USELOGSGLS){
	exp(@result$coefs)
}else{
	@result$coefs
}
@iconv <- @result$iconv
@iter <- @result$iter

if (@iconv == 0){
	print(paste("WARNING: failed to converge in",@maxit,"iterations"))
}
@eigs <- releigen(solve(@s), dmat(1/@psihat))
@J <- run(@m)
@vals <- 1/reverse(@eigs$values) - 1
@lhat <- if (@vals[@m] < 0){
	NULL
}else{
	@eigs$vectors[,@p + 1 - @J] * sqrt(@vals[@J]')
}
delete(@eigs,@J)

if (!@silent){
	if (isnull(@lhat)){
		print(paste("WARNING: eigenvalue",@m,"is negative"))
	}
	if (min(@psihat) < 0){
		print("WARNING: psi not acceptable since min(psihat) < 0")
	} elseif (min(eigenvals(@s - dmat(@psihat))) < 0) {
		print(paste("WARNING:",@what,\
			"- dmat(psihat) not non-negative definite"),macroname:T)
	}
	if (@iconv == 4){
		print("WARNING: could not reduce RSS on iteration",@iter)
	}elseif (@iconv <= 0){
		print(paste("WARNING: no convergence in",@iter,"iterations"))
	}
}
if (!@quiet && !@silent){
	if(@iconv > 0 && @iconv < 4){
		print(paste("Convergence in",@iter,"iterations by criterion",@iconv))
	}	
	print(psihat:@psihat,loadings:@lhat,criterion:@result$rss,vals:@vals)
}

delete(@quiet,@silent,@p,@m,@psi0,__USELOGSGLS)
@result <- structure(psihat:delete(@psihat,return:T),\
	loadings:delete(@lhat,return:T),criterion:delete(@result$rss,return:T),\
	vals:delete(@vals,return:T),status:delete(@iconv,return:T),\
	iter:delete(@iter,return:T))
delete(@result,return:T,invisible:T)
%glsfactor%

====> glsresids <====
glsresids   MACRO    DOLLARS
) Macro to be used with levmar to do GLS factor extraction
) This is very much a work in progress
) Usage:
)    stuff <- levmar(if(__USELOGSGLS){log(psi)}else{psi},\
)    vector(r),vector(r),m,resid:glsresids), m = number of factors
) glsresids itself returns vector(Ip - solve(r,sigmahat))
) It saves its most recent argument as _PSIGLS
)
) You need to set __USELOGSGLS to T or F before calling levmar
) and use log(psi) as an argument to levmar when __USELOGSGLS is True
) This is done automatically by macro glsfactor
)
) The result from levmar is the structure returned by levmar.
)   stuff$coefs is psihat (or log(psihat))
)   stuff$rss is the sum(vector(Ip - solve(r,sigmahat))^2)
)   stuff$iconv is convergence status -- 0 means not converged
)   stuff$iter is the number of iterations
) Version 991115
# $S(psi,vector(r),vector(r),m)
 # residualstr((Ip - solve(s,sigmahat))^2)
_PSIGLS <- if(__USELOGSGLS){
	exp($1)
}else{
	$1
}
@r <- $2 #x
@p <- sqrt(length(@r))
@r <- matrix(@r,@p)
@m <- $4
@eigs <- releigen(solve(@r), dmat(1/_PSIGLS))
@J <- run(@p, @p - @m + 1)
@tmpl <- @eigs$vectors[,@J]
@vals <- 1/@eigs$values[@J] - 1
@sighat <- (@tmpl * @vals') %C% @tmpl + dmat(_PSIGLS)
@resids <- sqrt(abs(dmat(@p,1) - solve(@r, @sighat)))
@resids <- @resids * @resids'
delete(@m,@J,@eigs,@tmpl,@vals,@p,@r,@sighat)
vector(delete(@resids, return:T))
%glsresids%

====> goodfit <====
goodfit             MACRO DOLLARS
) Usage: goodfit(s,lhat,psihat), s a covar or correl matrix
) Compute three goodness of fit criteria for factor analysis
) solution, the maximum absolute deviation, the sum of squared
) deviations and the sum of absolution deviations of the fitted
) covariance matrix lhat %*% lhat' + dmat(psihat) from the actual
) covariance matrix s
))Version 991114
# usage: $S(s,lhat,psihat), s a covar or correl matrix
@lhat <- argvalue($2,"lhat","real matrix nonmissing")
@dev <- vector(argvalue($1,"s","real square") -\
	(@lhat %C% @lhat)  -\
	dmat(argvalue($3[[1]],"psi","nonmissing real vector")))

@tmp <- vector(max(abs(@dev)),sum(@dev^2), sum(abs(@dev)),\
	labels:vector("max(|dev|)","sum(dev^2)","sum(|dev|)"))
delete(@dev)
delete(@tmp,return:T)
%goodfit%

====> stepuls <====
stepuls            MACRO DOLLARS
) usage: stepuls(s, psi, m [,print:T])
) One step of crude iteration for ULS factor analysis
) s = covariance or correlation matrix
) psi = vector of current values for uniquenesses or a structure
)       whose first component is such a vector
) m = number of factors
) If the optional 4th argument is T, the new values will be printed
) The output is a structure with components 'psi' and 'loadings'
) Version of 991114
# usage: $S(r,psi,m [,print:T]), s=covar or corr, psi=current value
# m = order = number of factors
if ($v != 3 || $k > 1){
	error("usage: $S(s, psi, m [print:T])")
}
@s <- matrix(argvalue($1,"s","real square"))
@psi <- argvalue($2[[1]],"psi","positive vector")
@M <- argvalue($3,"number of factors","positive integer scalar")
@print <- keyvalue($K,"print","TF",default:F)

@eigs <- eigen(@s - dmat(@psi))
@vals <- @eigs$values[run(@M)]
if(min(@vals)<= 0){
	error("s - psi not positive definite")
}
@l <- sqrt(@vals)' * @eigs$vectors[,run(@M)]
delete(@vals)
@psi <- diag(delete(@s,return:T)) - vector(sum((@l * @l)'))
if(min(@psi) <= 0){
	error("min(psi) <= 0")
}

CRITERION <- sum(@eigs$values[-run(@M)]^2)
delete(@eigs)
PSI <- delete(@psi, return:T)
LOADINGS <- delete(@l, return:T)
if(delete(@print,return:T)){ # print goodness of fit & psi
	print(psi:PSI ,criterion:CRITERION)
}
structure(psi:PSI, loadings:LOADINGS, crit:CRITERION)
%stepuls%

====> stepml <====
stepml             MACRO DOLLARS
) usage: stepml(s,psi,m [, print:T])
) One step of crude iteration for MLE factor analysis
) s = covariance or correlation matrix
) psi = vector of current values for uniquenesses or structure
)       whose first component is such a vector
) m = number of factors
) If the optional 4th argument is T, the new values will be printed
) The output is a structure with components 'psi' and 'loadings'
)) Version 991114
# usage: $S(s,psi,m [,print:T]), s=covar or corr, psi=current value
# m=order, print=T or F
if ($v != 3 || $k > 1){
	error("usage: $S(s, psi, m [print:T])")
}
@s <- matrix(argvalue($1,"argument 1","real square"))
@M <- argvalue($3,"number of factors","positive integer scalar")
@psi <- argvalue($2[[1]],"psi","positive vector")
@print <- keyvalue($K,"print","TF",default:F)

@eigs <- releigen(@s, dmat(@psi))
@vals <- @eigs$values[run(@M)] - 1
if(min(@vals)<= 0){
	error("s - psi not positive definite")
}
@l <- @psi * @eigs$vectors[,run(@M)] * sqrt(@vals)'
delete(@vals)

@psi <- diag(@s) - vector(sum((@l*@l)'))
if(min(@psi) <= 0){
	error("non-positive psi; cannot continue")
}
if (haslabels(@s)){
	@psi<-vector(@psi,labels:getlabels(@s,1))
}
CRITERION <- sum((@eigs$values - log(@eigs$values)-1)[-run(@M)])
delete(@eigs, @s)
PSI <- delete(@psi, return:T)
LOADINGS <- delete(@l, return:T)
if(delete(@print, return:T)){ # print goodness of fit & psi
	print(psi:PSI ,criterion:CRITERION)
}
structure(psi:PSI, loadings:LOADINGS, crit:CRITERION)
%stepml%

====> stepgls <====
stepgls            MACRO DOLLARS
) usage: stepgls(s, psi, m [, print:T])
) One step of crude iteration for GLS factor analysis
) s = covariance or correlation matrix
) psi = vector of current values for uniquenesses or a structure
)       whose first component is such a vector
) m = number of factors
) If the optional 4th argument is T, the new values will be printed
) The output is a structure with components 'psi' and 'loadings'
))Version 991114
# usage: $S(r, psi, m [,print:T]), s=covar or corr, psi=current value
# m = order = number of factors
if ($v != 3 || $k > 1){
	error("usage: $S(s, psi, m [print:T])")
}
@s <- matrix(argvalue($1,"argument 1","real square"))
@psi <- argvalue($2[[1]],"psi","positive vector")
@M <- argvalue($3,"number of factors","positive integer scalar")
@print <- keyvalue($K,"print","TF",default:F)

@eigs <- releigen(@s, dmat(@psi))
@vals <- @eigs$values[run(@M)] - 1
if(min(@vals)<= 0){
	error("s - psi not positive definite")
}
@l <- @psi * @eigs$vectors[,run(@M)] * sqrt(@vals)'

@sinv <- solve(@s)
@psi <- vector(solve(@sinv^2,\
	diag(@sinv) - vector(sum((@l %c% @sinv)^2))))
if(min(@psi) <= 0){
	error("non-positive psi; cannot continue")
}
CRITERION <- sum((1-1/@eigs$values[-run(@M)])^2)
delete(@eigs, @sinv, @s)
PSI <- delete(@psi, return:T)
LOADINGS <- delete(@l, return:T)
if(delete(@print,return:T)){ # print goodness of fit & psi
	print(psi:PSI ,criterion:CRITERION)
}
structure(psi:PSI, loadings:LOADINGS, crit:CRITERION)
%stepgls%