#
#   Copyright 2007-2021 by the individuals mentioned in the source code history
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#        http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.


library(OpenMx)
library(MASS)
library(testthat)
context("AlgebraCompute")

"%&%" <- OpenMx::"%&%"  # ensure we don't use the %&% from Matrix

A <- mxMatrix(values = runif(25), nrow = 5, ncol = 5, name = 'A')
B <- mxMatrix(values = runif(25), nrow = 5, ncol = 5, name = 'B')
C <- mxMatrix(values = runif(30), nrow = 6, ncol = 5, name = 'C')
D <- mxMatrix(values = 1:10, nrow = 2, ncol = 5, name = 'D')
E <- mxMatrix(values = 1:10, nrow = 5, ncol = 2, name = 'E')
G <- mxMatrix(values = c(0, -1, 1, -1), nrow=2, ncol=2, name='G')
Gprime <- mxMatrix(values = c(0, -1, 1, -1), nrow=2, ncol=2,
		   labels=paste0('D[1,',1:4,']'),
		   name='Gprime')
H <- mxMatrix(values = c(-9.0912244, 11.1561436, 12.9500773, -0.8400564, -5.2556798, -1.3265166, -4.8423977,
												 -14.4462163, -15.3022383, 17.1928667, -2.4835399, -18.1943775, 2.8587299, 18.6432244,
												 -13.5789122, -2.2622090), nrow=4, ncol=4, name='H')
# For Mnor
I <- mxMatrix("Iden", nrow = 5, ncol = 5, name="I")
L <- mxMatrix("Full", nrow=1, ncol=3, values=c(0,0,-Inf), free=F, name="L")
M <- mxMatrix("Full", nrow=1, ncol=3, values=c(0,0,0), free=F, name="M")
N <- mxMatrix("Full", nrow=1, ncol=3, values=c(0,1,0), free=F, name="N")
O <- mxMatrix("Full", nrow=3, ncol=3, values=1:9, free=F, name="O")
U <- mxMatrix("Full", nrow=1, ncol=3, values=c(Inf,Inf,Inf), free=F, name="U")
V <- mxMatrix("Stand", nrow=3, ncol=3, values=c(.5, .5, .5), free=F, name="V")
Z <- mxMatrix("Full", nrow=3, ncol=3, name='Z', values = c(6, 3, -1, 0, 8, 5, 0, 3, -1111))
EPSILON <- mxMatrix(values = 0.04*1:25, nrow = 5, ncol = 5, name = "EPSILON")
dimnames(A) <- list(letters[1:5], letters[22:26])
dimnames(B) <- dimnames(A)
dimnames(C) <- list(letters[1:6], letters[22:26])

model <- mxModel(A, B, C, D, E, G, Gprime, H, I, V, M, N, O, L, U, Z, EPSILON, name = 'model')
zeta <- 'z'
alpha <- 'a'
two <- 2

omxCheckTrue(all(mxEval(Gprime, model) == mxEval(G, model)))

# Insert passing tests
model <- mxModel(model, mxAlgebra(A * (B + A), name = 'parens'))
model <- mxModel(model, mxAlgebraFromString("A * (B + A)", name = 'parens2'))
model <- mxModel(model, mxAlgebra(A, name = 'test0'))
model <- mxModel(model, mxAlgebra(solve(A), name = 'test1'))
model <- mxModel(model, mxAlgebra(mpinv(C), name = 'test1a'))
model <- mxModel(model, mxAlgebra(t(A), name = 'test2'))
model <- mxModel(model, mxAlgebra(A ^ B, name = 'test3'))
model <- mxModel(model, mxAlgebra(t(A) ^ B, name = 'test3a'))
model <- mxModel(model, mxAlgebra(A ^ t(B), name = 'test3b'))
model <- mxModel(model, mxAlgebra(A %*% B, name = 'test4'))
model <- mxModel(model, mxAlgebra(t(A) %*% B, name = 'test4a'))
model <- mxModel(model, mxAlgebra(A %*% t(B), name = 'test4b'))
model <- mxModel(model, mxAlgebra(A * B, name = 'test5'))
model <- mxModel(model, mxAlgebra(t(A) * B, name = 'test5a'))
model <- mxModel(model, mxAlgebra(A * t(B), name = 'test5b'))
model <- mxModel(model, mxAlgebra(t(A) * t(B), name = 'test5c'))
model <- mxModel(model, mxAlgebra(A %x% B, name = 'test6'))
model <- mxModel(model, mxAlgebra(t(A) %x% B, name = 'test6a'))
model <- mxModel(model, mxAlgebra(A %x% t(B), name = 'test6b'))
model <- mxModel(model, mxAlgebra(A %^% B, name = 'test38'))
model <- mxModel(model, mxAlgebra(t(A) %^% B, name = 'test38a'))
model <- mxModel(model, mxAlgebra(A %^% t(B), name = 'test38b'))
model <- mxModel(model, mxAlgebra(A %&% B, name = 'test7'))
model <- mxModel(model, mxAlgebra(t(A) %&% B, name = 'test7a'))
model <- mxModel(model, mxAlgebra(A %&% t(B), name = 'test7b'))
model <- mxModel(model, mxAlgebra(C %&% B, name = 'test7c'))
model <- mxModel(model, mxAlgebra(C %&% t(B), name = 'test7d'))
model <- mxModel(model, mxAlgebra(A / B, name = 'test8'))
model <- mxModel(model, mxAlgebra(t(A) / B, name = 'test8a'))
model <- mxModel(model, mxAlgebra(A / t(B), name = 'test8b'))
model <- mxModel(model, mxAlgebra(t(A) / t(B), name = 'test8c'))
model <- mxModel(model, mxAlgebra(A + B, name = 'test9'))
model <- mxModel(model, mxAlgebra(t(A) + B, name = 'test9a'))
model <- mxModel(model, mxAlgebra(A + t(B), name = 'test9b'))
model <- mxModel(model, mxAlgebra(t(A) + t(B), name = 'test9c'))
model <- mxModel(model, mxAlgebra(A - B, name = 'test10'))
model <- mxModel(model, mxAlgebra(t(A) - B, name = 'test10a'))
model <- mxModel(model, mxAlgebra(A - t(B), name = 'test10b'))
model <- mxModel(model, mxAlgebra(t(A) - t(B), name = 'test10c'))
model <- mxModel(model, mxAlgebra(-A, name = 'test11'))
model <- mxModel(model, mxAlgebra(-t(A), name = 'test11a'))
model <- mxModel(model, mxAlgebra(t(-A), name = 'test11b'))
model <- mxModel(model, mxAlgebra(cbind(A,B), name = 'test12'))
model <- mxModel(model, mxAlgebra(cbind(t(A),B), name = 'test12a'))
model <- mxModel(model, mxAlgebra(cbind(A,t(B)), name = 'test12b'))
model <- mxModel(model, mxAlgebra(rbind(A,B), name = 'test13'))
model <- mxModel(model, mxAlgebra(rbind(t(A),B), name = 'test13a'))
model <- mxModel(model, mxAlgebra(rbind(A,t(B)), name = 'test13b'))
model <- mxModel(model, mxAlgebra(det(A), name = 'test14'))
model <- mxModel(model, mxAlgebra(det(t(A)), name = 'test14a'))
model <- mxModel(model, mxAlgebra(tr(A), name = 'test15'))
model <- mxModel(model, mxAlgebra(sum(A,B), name = 'test16'))
model <- mxModel(model, mxAlgebra(prod(A,B), name = 'test17'))
model <- mxModel(model, mxAlgebra(max(A,B), name = 'test18'))
model <- mxModel(model, mxAlgebra(min(A,B), name = 'test19'))
model <- mxModel(model, mxAlgebra(abs(A), name = 'test20'))
model <- mxModel(model, mxAlgebra(cos(A), name = 'test21'))
model <- mxModel(model, mxAlgebra(cosh(A), name = 'test22'))
model <- mxModel(model, mxAlgebra(sin(A), name = 'test23'))
model <- mxModel(model, mxAlgebra(sinh(A), name = 'test24'))
model <- mxModel(model, mxAlgebra(tan(A), name = 'test25'))
model <- mxModel(model, mxAlgebra(tanh(A), name = 'test26'))
model <- mxModel(model, mxAlgebra(exp(A), name = 'test27'))
model <- mxModel(model, mxAlgebra(log(A), name = 'test28'))
model <- mxModel(model, mxAlgebra(sqrt(A), name = 'test29'))
model <- mxModel(model, mxAlgebra(A[1,1], name = 'test30a'))
model <- mxModel(model, mxAlgebra(A[1,0], name = 'test30b'))
model <- mxModel(model, mxAlgebra(A[0,1], name = 'test30c'))
model <- mxModel(model, mxAlgebra(A[0,0], name = 'test30d'))
model <- mxModel(model, mxAlgebra(A[,], name = 'test30e'))
model <- mxModel(model, mxAlgebra(A[1,], name = 'test30f'))
model <- mxModel(model, mxAlgebra(A[,1], name = 'test30g'))
model <- mxModel(model, mxAlgebra(A[0,], name = 'test30h'))
model <- mxModel(model, mxAlgebra(A[,0], name = 'test30i'))
model <- mxModel(model, mxAlgebra(A[,1] + B[,2], name = 'test30j'))
model <- mxModel(model, mxAlgebra(A[1,] + B[2,], name = 'test30k'))
model <- mxModel(model, mxAlgebra(A['a','v'], name = 'test30l'))
model <- mxModel(model, mxAlgebra(A['d','v'], name = 'test30m'))
model <- mxModel(model, mxAlgebra(A['e',], name = 'test30n'))
model <- mxModel(model, mxAlgebra(A[,'z'], name = 'test30o'))
model <- mxModel(model, mxAlgebra(A[,zeta], name = 'test30p'))
model <- mxModel(model, mxAlgebra(A[alpha,], name = 'test30q'))
model <- mxModel(model, mxAlgebra(A[two,two], name = 'test30r'))
model <- mxModel(model, mxAlgebra(A[1:3,], name = 'test30s'))
model <- mxModel(model, mxAlgebra(A[,2:4], name = 'test30t'))
model <- mxModel(model, mxAlgebra(A[cbind(1,1), cbind(2,2)], name = 'test30u'))
model <- mxModel(model, mxAlgebra(A[-1:-3, -2], name = 'test30v'))
model <- mxModel(model, mxAlgebra(A[-2:-4,3], name = 'test30w'))
model <- mxModel(model, mxAlgebra(A[3, -2:-4], name = 'test30x'))
model <- mxModel(model, mxAlgebra(A[cbind('a','b','c'), cbind('v','z')], name = 'test30y'))
model <- mxModel(model, mxAlgebra(vech(D), name = 'test31a'))
model <- mxModel(model, mxAlgebra(vech(E), name = 'test31b'))
model <- mxModel(model, mxAlgebra(vechs(D), name = 'test32a'))
model <- mxModel(model, mxAlgebra(vechs(E), name = 'test32b'))
model <- mxModel(model, mxAlgebra(diag2vec(B), name = 'test33a'))
model <- mxModel(model, mxAlgebra(diag2vec(C), name = 'test33b'))
model <- mxModel(model, mxAlgebra(vec2diag(B[1,]), name = 'test34a'))
model <- mxModel(model, mxAlgebra(vec2diag(t(B[1,])), name = 'test34b'))
model <- mxModel(model, mxAlgebra(omxMnor(V, M, L, U), name= 'test35a'))
model <- mxModel(model, mxAlgebra(3 : 5.2, name = 'test37a'))
model <- mxModel(model, mxAlgebra(5.2 : 3, name = 'test37b'))
model <- mxModel(model, mxAlgebra(0 : 0, name = 'test37c'))
model <- mxModel(model, mxAlgebra(-3 : -5.2, name = 'test37d'))
model <- mxModel(model, mxAlgebra(-5.2 : -3, name = 'test37e'))
model <- mxModel(model, mxAlgebra(rvectorize(A), name = 'test39a'))
model <- mxModel(model, mxAlgebra(rvectorize(t(A)), name = 'test39b'))
model <- mxModel(model, mxAlgebra(cvectorize(A), name = 'test39c'))
model <- mxModel(model, mxAlgebra(cvectorize(t(A)), name = 'test39d'))
model <- mxModel(model, mxAlgebra(cvectorize(rvectorize(A)), name = 'test39e'))
model <- mxModel(model, mxAlgebra(cvectorize(rvectorize(t(A))), name = 'test39f'))
model <- mxModel(model, mxAlgebra(rvectorize(cvectorize(A)), name = 'test39g'))
model <- mxModel(model, mxAlgebra(rvectorize(cvectorize(t(A))), name = 'test39h'))
model <- mxModel(model, mxAlgebra(eigenvec(A), name='test40'))
model <- mxModel(model, mxAlgebra(eigenvec(G), name='test40a'))
model <- mxModel(model, mxAlgebra(eigenvec(H), name='test40b'))
model <- mxModel(model, mxAlgebra(eigenval(A), name='test41'))
model <- mxModel(model, mxAlgebra(eigenval(G), name='test41a'))
model <- mxModel(model, mxAlgebra(eigenval(H), name='test41b'))
model <- mxModel(model, mxAlgebra(ieigenvec(A), name='test42'))
model <- mxModel(model, mxAlgebra(ieigenvec(G), name='test42a'))
model <- mxModel(model, mxAlgebra(ieigenvec(H), name='test42b'))
model <- mxModel(model, mxAlgebra(ieigenval(A), name='test43'))
model <- mxModel(model, mxAlgebra(ieigenval(G), name='test43a'))
model <- mxModel(model, mxAlgebra(ieigenval(H), name='test43b'))
model <- mxModel(model, mxAlgebra(omxNot(Z), name ='test44a'))
model <- mxModel(model, mxAlgebra(omxNot(omxNot(Z)), name = 'test44b'))
model <- mxModel(model, mxAlgebra(omxNot(t(Z)), name = 'test44c'))
model <- mxModel(model, mxAlgebra(omxSelectRows(O, N), name='test45'))
model <- mxModel(model, mxAlgebra(omxSelectCols(O, N), name='test46'))
model <- mxModel(model, mxAlgebra(omxSelectRowsAndCols(O, N), name='test47'))
model <- mxModel(model, mxAlgebra(omxSelectRows(O,omxNot(N)), name='test45b'))
model <- mxModel(model, mxAlgebra(omxSelectCols(O,omxNot(N)), name='test46b'))
model <- mxModel(model, mxAlgebra(omxSelectRowsAndCols(O,omxNot(N)), name='test47b'))
model <- mxModel(model, mxAlgebra(omxSelectRows(O,N%x%5), name='test45c'))
model <- mxModel(model, mxAlgebra(omxSelectCols(O,N%x%5), name='test46c'))
model <- mxModel(model, mxAlgebra(omxSelectRowsAndCols(O,N%x%5), name='test47c'))
model <- mxModel(model, mxAlgebra(omxSelectRows(O,t(N)), name='test45d'))
model <- mxModel(model, mxAlgebra(omxSelectCols(O,t(N)), name='test46d'))
model <- mxModel(model, mxAlgebra(omxSelectRowsAndCols(O,t(N)), name='test47d'))
model <- mxModel(model, mxAlgebra(omxSelectRows(O,omxNot(t(N))), name='test45e'))
model <- mxModel(model, mxAlgebra(omxSelectCols(O,omxNot(t(N))), name='test46e'))
model <- mxModel(model, mxAlgebra(omxSelectRowsAndCols(O,omxNot(t(N))), name='test47e'))
model <- mxModel(model, mxAlgebra(omxSelectRows(O,t(N)%x%5), name='test45f'))
model <- mxModel(model, mxAlgebra(omxSelectCols(O,t(N)%x%5), name='test46f'))
model <- mxModel(model, mxAlgebra(omxSelectRowsAndCols(O,t(N)%x%5), name='test47f'))
model <- mxModel(model, mxAlgebra(omxSelectRows(O[,1:2],N%x%5),    name='test45g'))
model <- mxModel(model, mxAlgebra(omxSelectCols(O[1:2,],N%x%5),    name='test46g'))
model <- mxModel(model, mxAlgebra(omxSelectRows(O[,1:2],t(N)%x%5), name='test45h'))
model <- mxModel(model, mxAlgebra(omxSelectCols(O[1:2,],t(N)%x%5), name='test46h'))
model <- mxModel(model, mxAlgebra(mean(A), name='test48a'))
model <- mxModel(model, mxAlgebra(mean(B), name='test48b'))
model <- mxModel(model, mxAlgebra(mean(C), name='test48c'))
model <- mxModel(model, mxAlgebra(mean(t(C)), name='test48d'))
model <- mxModel(model, mxAlgebra(mean(-C), name='test48e'))
model <- mxModel(model, mxAlgebra(omxGreaterThan(A,B), name ='test49a'))
model <- mxModel(model, mxAlgebra(omxGreaterThan(B,A), name ='test49b'))
model <- mxModel(model, mxAlgebra(omxGreaterThan(A,A), name ='test49c'))
model <- mxModel(model, mxAlgebra(omxGreaterThan(B,B), name ='test49d'))
model <- mxModel(model, mxAlgebra(omxLessThan(A,B), name ='test50a'))
model <- mxModel(model, mxAlgebra(omxLessThan(B,A), name ='test50b'))
model <- mxModel(model, mxAlgebra(omxLessThan(A,A), name ='test50c'))
model <- mxModel(model, mxAlgebra(omxLessThan(B,B), name ='test50d'))
model <- mxModel(model, mxAlgebra(omxAnd(omxNot(A),B), name ='test51a'))
model <- mxModel(model, mxAlgebra(omxAnd(B,A), name ='test51b'))
model <- mxModel(model, mxAlgebra(omxAnd(omxNot(A),A), name ='test51c'))
model <- mxModel(model, mxAlgebra(omxAnd(A,omxNot(A)), name ='test51d'))
model <- mxModel(model, mxAlgebra(omxOr(omxNot(A),B), name ='test52a'))
model <- mxModel(model, mxAlgebra(omxOr(B,A), name ='test52b'))
model <- mxModel(model, mxAlgebra(omxOr(omxNot(A),omxNot(A)), name ='test52c'))
model <- mxModel(model, mxAlgebra(omxOr(A,A), name ='test52d'))
model <- mxModel(model, mxAlgebra(omxApproxEquals(A,B, EPSILON), name ='test53a'))
model <- mxModel(model, mxAlgebra(omxApproxEquals(B,A, EPSILON), name ='test53b'))
model <- mxModel(model, mxAlgebra(omxApproxEquals(A,A, EPSILON), name ='test53c'))
model <- mxModel(model, mxAlgebra(omxApproxEquals(B,B, EPSILON), name ='test53d'))
model <- mxModel(model, mxAlgebra(omxApproxEquals(A,B,t(EPSILON)), name ='test53e'))
model <- mxModel(model, mxAlgebra(omxExponential(A), name = 'test55a'))
model <- mxModel(model, mxAlgebra(omxExponential(t(A)), name = 'test55b'))
model <- mxModel(model, mxAlgebra(omxExponential(B), name = 'test55c'))
model <- mxModel(model, mxAlgebra(omxExponential(t(B)), name = 'test55d'))
model <- mxModel(model, mxAlgebra(omxExponential(I), name = 'test55e'))
model <- mxModel(model, mxAlgebra(omxExponential(t(I)), name = 'test55f'))
model <- mxModel(model, mxAlgebra(chol(V), name = 'test56'))
model <- mxModel(model, mxAlgebra(cov2cor(A), name = 'test57a'))
model <- mxModel(model, mxAlgebra(cov2cor(B), name = 'test57b'))
model <- mxModel(model, mxAlgebra(vech2full( vech( t(A) %*% A ) ), name = 'test58a'))
model <- mxModel(model, mxAlgebra(vechs2full( vechs( cov2cor(t(A) %*% A) ) ), name = 'test58b'))
model <- mxModel(model, mxAlgebra(omxSelectRows(t(O), N), name='test65'))
model <- mxModel(model, mxAlgebra(omxSelectCols(t(O), N), name='test66'))
model <- mxModel(model, mxAlgebra(omxSelectRowsAndCols(t(O), N), name='test67'))
model <- mxModel(model, mxAlgebra(omxSelectRows(t(O),omxNot(N)), name='test65b'))
model <- mxModel(model, mxAlgebra(omxSelectCols(t(O),omxNot(N)), name='test66b'))
model <- mxModel(model, mxAlgebra(omxSelectRowsAndCols(t(O),omxNot(N)), name='test67b'))
model <- mxModel(model, mxAlgebra(omxSelectRows(t(O),N%x%5), name='test65c'))
model <- mxModel(model, mxAlgebra(omxSelectCols(t(O),N%x%5), name='test66c'))
model <- mxModel(model, mxAlgebra(omxSelectRowsAndCols(t(O),N%x%5), name='test67c'))
model <- mxModel(model, mxAlgebra(omxSelectRows(t(O),t(N)), name='test65d'))
model <- mxModel(model, mxAlgebra(omxSelectCols(t(O),t(N)), name='test66d'))
model <- mxModel(model, mxAlgebra(omxSelectRowsAndCols(t(O),t(N)), name='test67d'))
model <- mxModel(model, mxAlgebra(omxSelectRows(t(O),omxNot(t(N))), name='test65e'))
model <- mxModel(model, mxAlgebra(omxSelectCols(t(O),omxNot(t(N))), name='test66e'))
model <- mxModel(model, mxAlgebra(omxSelectRowsAndCols(t(O),omxNot(t(N))), name='test67e'))
model <- mxModel(model, mxAlgebra(omxSelectRows(t(O),t(N)%x%5), name='test65f'))
model <- mxModel(model, mxAlgebra(omxSelectCols(t(O),t(N)%x%5), name='test66f'))
model <- mxModel(model, mxAlgebra(omxSelectRowsAndCols(t(O),t(N)%x%5), name='test67f'))
model <- mxModel(model, mxAlgebra(logm(A), name = 'test68a'))
model <- mxModel(model, mxAlgebra(logm(B), name = 'test68b'))
model <- mxModel(model, mxAlgebra(logm(I), name = 'test68c'))
model <- mxModel(model, mxAlgebra(asin(A),name="test69a"))
model <- mxModel(model, mxAlgebra(acos(A),name="test69b"))
model <- mxModel(model, mxAlgebra(atan(A),name="test69c"))
model <- mxModel(model, mxAlgebra(lgamma(A),name="test70"))
model <- mxModel(model, mxAlgebra(p2z(A),name="test71a"))
model <- mxModel(model, mxAlgebra(asinh(A),name="test72a"))
model <- mxModel(model, mxAlgebra(acosh(D),name="test72b"))
model <- mxModel(model, mxAlgebra(atanh(A),name="test72c"))
model <- mxModel(model, mxAlgebra(lgamma1p(A),name="test73"))
model <- mxModel(model, mxAlgebra(logp2z(log(A)),name="test71b"))
model <- mxModel(model, mxAlgebra(dbeta(A,A,A,A,A),name="test74a"))
model <- mxModel(model, mxAlgebra(pbeta(A,A,A,A,A,A),name="test75a"))
model <- mxModel(model, mxAlgebra(dbeta(A,1,1,-1,0),name="test74b"))
model <- mxModel(model, mxAlgebra(pbeta(A,1,1,-1,1,0),name="test75b"))
model <- mxModel(model, mxAlgebra(dbeta(C,D,D,-1,0),name="test74c"))
model <- mxModel(model, mxAlgebra(pbeta(C,D,D,-1,1,0),name="test75c"))
model <- mxModel(model, mxAlgebra(besselI(A,A,A),name="test76a"))
model <- mxModel(model, mxAlgebra(besselI(C,D,0),name="test76b"))
model <- mxModel(model, mxAlgebra(besselJ(A,A),name="test77a"))
model <- mxModel(model, mxAlgebra(besselJ(C,D),name="test77b"))
model <- mxModel(model, mxAlgebra(besselK(A,A,A),name="test78a"))
model <- mxModel(model, mxAlgebra(besselK(C,D,0),name="test78b"))
model <- mxModel(model, mxAlgebra(besselY(A,A),name="test79a"))
model <- mxModel(model, mxAlgebra(besselY(C,D),name="test79b"))
model <- mxModel(model, mxAlgebra(dpois(D,D,0),name="test80a"))
model <- mxModel(model, mxAlgebra(dpois(D,D,D),name="test80b"))
model <- mxModel(model, mxAlgebra(ppois(D,D,D,D),name="test81a"))
model <- mxModel(model, mxAlgebra(ppois(D,D,1,0),name="test81b"))
model <- mxModel(model, mxAlgebra(ppois(D,D,0,0),name="test81c"))
model <- mxModel(model, mxAlgebra(omxDnbinom(D,D,C,-1,0),name="test82a"))
model <- mxModel(model, mxAlgebra(omxDnbinom(D,D,C,-1,C),name="test82b"))
model <- mxModel(model, mxAlgebra(omxPnbinom(D,D,C,-1,C,0),name="test83a"))
model <- mxModel(model, mxAlgebra(omxPnbinom(D,D,C,-1,1,0),name="test83b"))
model <- mxModel(model, mxAlgebra(omxPnbinom(D,D,C,-1,0,1),name="test83c"))
model <- mxModel(model, mxAlgebra(omxDnbinom(D,D,-1,C,0),name="test84a"))
model <- mxModel(model, mxAlgebra(omxDnbinom(D,D,-1,C,C),name="test84b"))
model <- mxModel(model, mxAlgebra(omxPnbinom(D,D,-1,C,C,0),name="test85a"))
model <- mxModel(model, mxAlgebra(omxPnbinom(D,D,-1,C,1,0),name="test85b"))
model <- mxModel(model, mxAlgebra(omxPnbinom(D,D,-1,C,0,1),name="test85c"))
model <- mxModel(model, mxAlgebra(dchisq(A,A,A,A),name="test86a"))
model <- mxModel(model, mxAlgebra(pchisq(A,A,A,A,A),name="test86b"))
model <- mxModel(model, mxAlgebra(dchisq(D,1,-1,0),name="test87a"))
model <- mxModel(model, mxAlgebra(pchisq(D,1,-1,0,0),name="test87b"))
model <- mxModel(model, mxAlgebra(dbinom(D,D,A,0),name="test88a"))
model <- mxModel(model, mxAlgebra(dbinom(D,D,A,1),name="test88b"))
model <- mxModel(model, mxAlgebra(pbinom(D,D,A,0,0),name="test89a"))
model <- mxModel(model, mxAlgebra(pbinom(D,D,A,1,1),name="test89b"))
model <- mxModel(model, mxAlgebra(dcauchy(H,A,D,0),name="test90a"))
model <- mxModel(model, mxAlgebra(dcauchy(H,A,D,1),name="test90b"))
model <- mxModel(model, mxAlgebra(pcauchy(H,A,D,0,1),name="test91a"))
model <- mxModel(model, mxAlgebra(pcauchy(H,A,D,1,0),name="test91b"))
model <- mxModel(model, mxAlgebra(rowSums(C),name="test92"))
model <- mxModel(model, mxAlgebra(colSums(C),name="test93"))
model <- mxModel(model, mxAlgebra(1/(1+N), name='test94'))
modelOut <- mxRun(model)

# Check passing tests
omxCheckCloseEnough(modelOut[['test0']]$result, A$values, 0.001)
omxCheckCloseEnough(modelOut[['parens']]$result, A$values * (B$values + A$values), 0.001)
omxCheckCloseEnough(modelOut[['parens2']]$result, A$values * (B$values + A$values), 0.001)
omxCheckCloseEnough(modelOut[['test1']]$result, solve(A$values), 0.001)
omxCheckCloseEnough(modelOut[['test1a']]$result, mpinv(C$values), 0.001)
omxCheckCloseEnough(modelOut[['test2']]$result, t(A$values), 0.001)
omxCheckCloseEnough(modelOut[['test3']]$result, A$values ^ B$values, 0.001)
omxCheckCloseEnough(modelOut[['test3a']]$result, t(A$values) ^ B$values, 0.001)
omxCheckCloseEnough(modelOut[['test3b']]$result, A$values ^ t(B$values), 0.001)
omxCheckCloseEnough(modelOut[['test4']]$result, A$values %*% B$values, 0.001)
omxCheckCloseEnough(modelOut[['test4a']]$result, t(A$values) %*% B$values, 0.001)
omxCheckCloseEnough(modelOut[['test4b']]$result, A$values %*% t(B$values), 0.001)
omxCheckCloseEnough(modelOut[['test5']]$result, A$values * B$values, 0.001)
omxCheckCloseEnough(modelOut[['test5a']]$result, t(A$values) * B$values, 0.001)
omxCheckCloseEnough(modelOut[['test5b']]$result, A$values * t(B$values), 0.001)
omxCheckCloseEnough(modelOut[['test5c']]$result, t(A$values) * t(B$values), 0.001)
omxCheckCloseEnough(modelOut[['test6']]$result, A$values %x% B$values, 0.001)
omxCheckCloseEnough(modelOut[['test6a']]$result, t(A$values) %x% B$values, 0.001)
omxCheckCloseEnough(modelOut[['test6b']]$result, A$values %x% t(B$values), 0.001)
omxCheckCloseEnough(modelOut[['test38']]$result, A$values %^% B$values, 0.001)
omxCheckCloseEnough(modelOut[['test38a']]$result, t(A$values) %^% B$values, 0.001)
omxCheckCloseEnough(modelOut[['test38b']]$result, A$values %^% t(B$values), 0.001)
omxCheckCloseEnough(modelOut[['test7']]$result, A$values %&% B$values, 0.001)
omxCheckCloseEnough(modelOut[['test7a']]$result, t(A$values) %&% B$values, 0.001)
omxCheckCloseEnough(modelOut[['test7b']]$result, A$values %&% t(B$values), 0.001)
omxCheckCloseEnough(modelOut[['test7c']]$result, C$values %&% B$values, 0.001)
omxCheckCloseEnough(modelOut[['test7d']]$result, C$values %&% t(B$values), 0.001)
omxCheckCloseEnough(modelOut[['test8']]$result, A$values / B$values, 0.001)
omxCheckCloseEnough(modelOut[['test8a']]$result, t(A$values) / B$values, 0.001)
omxCheckCloseEnough(modelOut[['test8b']]$result, A$values / t(B$values), 0.001)
omxCheckCloseEnough(modelOut[['test8c']]$result, t(A$values) / t(B$values), 0.001)
omxCheckCloseEnough(modelOut[['test9']]$result, A$values + B$values, 0.001)
omxCheckCloseEnough(modelOut[['test9a']]$result, t(A$values) + B$values, 0.001)
omxCheckCloseEnough(modelOut[['test9b']]$result, A$values + t(B$values), 0.001)
omxCheckCloseEnough(modelOut[['test9c']]$result, t(A$values) + t(B$values), 0.001)
omxCheckCloseEnough(modelOut[['test10']]$result, A$values - B$values, 0.001)
omxCheckCloseEnough(modelOut[['test10a']]$result, t(A$values) - B$values, 0.001)
omxCheckCloseEnough(modelOut[['test10b']]$result, A$values - t(B$values), 0.001)
omxCheckCloseEnough(modelOut[['test10c']]$result, t(A$values) - t(B$values), 0.001)
omxCheckCloseEnough(modelOut[['test11']]$result, -A$values, 0.001)
omxCheckCloseEnough(modelOut[['test11a']]$result, -t(A$values), 0.001)
omxCheckCloseEnough(modelOut[['test11b']]$result, t(-A$values), 0.001)
omxCheckCloseEnough(mxEval(test12, modelOut), cbind(A$values, B$values), 0.001)
omxCheckCloseEnough(mxEval(test12a, modelOut), cbind(t(A$values), B$values), 0.001)
omxCheckCloseEnough(mxEval(test12b, modelOut), cbind(A$values, t(B$values)), 0.001)
omxCheckCloseEnough(mxEval(test13, modelOut), rbind(A$values, B$values), 0.001)
omxCheckCloseEnough(mxEval(test13a, modelOut), rbind(t(A$values), B$values), 0.001)
omxCheckCloseEnough(mxEval(test13b, modelOut), rbind(A$values, t(B$values)), 0.001)
omxCheckCloseEnough(modelOut[['test14']]$result, det(A$values), 0.001)
omxCheckCloseEnough(modelOut[['test14a']]$result, det(t(A$values)), 0.001)
omxCheckCloseEnough(modelOut[['test15']]$result, tr(A$values), 0.001)
omxCheckCloseEnough(modelOut[['test16']]$result, sum(A$values, B$values), 0.001)
omxCheckCloseEnough(modelOut[['test17']]$result, prod(A$values, B$values), 0.001)
omxCheckCloseEnough(modelOut[['test18']]$result, max(A$values, B$values), 0.001)
omxCheckCloseEnough(modelOut[['test19']]$result, min(A$values, B$values), 0.001)
omxCheckCloseEnough(modelOut[['test20']]$result, abs(A$values), 0.001)
omxCheckCloseEnough(modelOut[['test21']]$result, cos(A$values), 0.001)
omxCheckCloseEnough(modelOut[['test22']]$result, cosh(A$values), 0.001)
omxCheckCloseEnough(modelOut[['test23']]$result, sin(A$values), 0.001)
omxCheckCloseEnough(modelOut[['test24']]$result, sinh(A$values), 0.001)
omxCheckCloseEnough(modelOut[['test25']]$result, tan(A$values), 0.001)
omxCheckCloseEnough(modelOut[['test26']]$result, tanh(A$values), 0.001)
omxCheckCloseEnough(modelOut[['test27']]$result, exp(A$values), 0.001)
omxCheckCloseEnough(modelOut[['test28']]$result, log(A$values), 0.001)
omxCheckCloseEnough(modelOut[['test29']]$result, sqrt(A$values), 0.001)
omxCheckCloseEnough(modelOut[['test30a']]$result, A$values[1,1], 0.001)
suppressWarnings(omxCheckCloseEnough(modelOut[['test30b']]$result, A$values[1,0], 0.001))
suppressWarnings(omxCheckCloseEnough(modelOut[['test30c']]$result, A$values[0,1], 0.001))
suppressWarnings(omxCheckCloseEnough(modelOut[['test30d']]$result, A$values[0,0], 0.001))
omxCheckCloseEnough(modelOut[['test30e']]$result, mxEval(A[,], model), 0.001)
omxCheckCloseEnough(modelOut[['test30f']]$result, mxEval(A[1,], model), 0.001)
omxCheckCloseEnough(modelOut[['test30g']]$result, mxEval(A[,1], model), 0.001)
omxCheckEquals(dim(modelOut[['test30h']]$result), c(0,0))
omxCheckEquals(dim(modelOut[['test30i']]$result), c(0,0))
omxCheckCloseEnough(modelOut[['test30j']]$result, mxEval(A[,1] + B[,2], model), 0.001)
omxCheckCloseEnough(modelOut[['test30k']]$result, mxEval(A[1,] + B[2,], model), 0.001)
omxCheckCloseEnough(modelOut[['test30l']]$result, A$values['a','v'], 0.001)
omxCheckCloseEnough(modelOut[['test30m']]$result, A$values['d','v'], 0.001)
omxCheckCloseEnough(modelOut[['test30n']]$result, mxEval(A['e',], model), 0.001)
omxCheckCloseEnough(modelOut[['test30o']]$result, mxEval(A[,'z'], model), 0.001)
omxCheckCloseEnough(modelOut[['test30p']]$result, mxEval(A[,zeta], model), 0.001)
omxCheckCloseEnough(modelOut[['test30q']]$result, mxEval(A[alpha,], model), 0.001)
omxCheckCloseEnough(modelOut[['test30r']]$result, mxEval(A[two,two], model), 0.001)
omxCheckCloseEnough(modelOut[['test30s']]$result, mxEval(A[1:3,], model), 0.001)
omxCheckCloseEnough(modelOut[['test30t']]$result, mxEval(A[,2:4], model), 0.001)
omxCheckCloseEnough(modelOut[['test30u']]$result, mxEval(A[cbind(1,1), cbind(2,2)], model), 0.001)
omxCheckCloseEnough(modelOut[['test30v']]$result, mxEval(A[-1:-3, -2], model), 0.001)
omxCheckCloseEnough(modelOut[['test30w']]$result, mxEval(A[-2:-4, 3], model), 0.001)
omxCheckCloseEnough(modelOut[['test30x']]$result, mxEval(A[3, -2:-4], model), 0.001)
omxCheckCloseEnough(modelOut[['test30y']]$result, mxEval(A[cbind('a','b','c'), cbind('v','z')], model), 0.001)
omxCheckCloseEnough(modelOut[['test31a']]$result, mxEval(as.matrix(vech(D)), model), 0.001)
omxCheckCloseEnough(modelOut[['test31b']]$result, mxEval(as.matrix(vech(E)), model), 0.001)
omxCheckCloseEnough(modelOut[['test32a']]$result, mxEval(as.matrix(vechs(D)), model), 0.001)
omxCheckCloseEnough(modelOut[['test32b']]$result, mxEval(as.matrix(vechs(E)), model), 0.001)
omxCheckCloseEnough(modelOut[['test33a']]$result, mxEval(diag2vec(B), model), 0.001)
omxCheckCloseEnough(modelOut[['test33b']]$result, mxEval(diag2vec(C), model), 0.001)
omxCheckCloseEnough(modelOut[['test34a']]$result, mxEval(vec2diag(B[1,]), model), 0.001)
omxCheckCloseEnough(modelOut[['test34b']]$result, mxEval(vec2diag(t(B[1,])), model), 0.001)
omxCheckCloseEnough(modelOut[['test35a']]$result, .33333, 0.001)
omxCheckCloseEnough(modelOut[['test37a']]$result, as.matrix(3:5.2), 0.001)
omxCheckCloseEnough(modelOut[['test37b']]$result, as.matrix(5.2:3), 0.001)
omxCheckCloseEnough(modelOut[['test37c']]$result, as.matrix(0:0), 0.001)
omxCheckCloseEnough(modelOut[['test37d']]$result, as.matrix(-3:-5.2), 0.001)
omxCheckCloseEnough(modelOut[['test37e']]$result, as.matrix(-5.2:-3), 0.001)
omxCheckCloseEnough(modelOut[['test39a']]$result, rvectorize(A$values), 0.001)
omxCheckCloseEnough(modelOut[['test39b']]$result, rvectorize(t(A$values)), 0.001)
omxCheckCloseEnough(modelOut[['test39c']]$result, cvectorize(A$values), 0.001)
omxCheckCloseEnough(modelOut[['test39d']]$result, cvectorize(t(A$values)), 0.001)
omxCheckCloseEnough(modelOut[['test39e']]$result, cvectorize(rvectorize(A$values)), 0.001)
omxCheckCloseEnough(modelOut[['test39f']]$result, cvectorize(rvectorize(t(A$values))), 0.001)
omxCheckCloseEnough(modelOut[['test39g']]$result, rvectorize(cvectorize(A$values)), 0.001)
omxCheckCloseEnough(modelOut[['test39h']]$result, rvectorize(cvectorize(t(A$values))), 0.001)
omxCheckCloseEnough(modelOut[["Gprime"]]$values, mxEval(Gprime, model, compute=T), 1e-4)

check_eigen_decomp <- function(mat, suffix) {
  algName <- paste0('test', 40:43, suffix)
  vec <- (modelOut[[ algName[1] ]]$result +
            1i * modelOut[[ algName[3] ]]$result)
  val <- (modelOut[[ algName[2] ]]$result +
            1i * modelOut[[ algName[4] ]]$result)
  z <- matrix(as.complex(0), nrow(mat), ncol(mat))
  expect_equivalent(vec %*% diag(c(val)) %*% solve(vec) - mat, z,
                    tolerance=1e-13, scale=1)
}

check_eigen_decomp(A$values, '')
check_eigen_decomp(G$values, 'a')
check_eigen_decomp(H$values, 'b')

omxCheckCloseEnough(modelOut[['test44a']]$result, mxEval(omxNot(Z), model), 0.001)
omxCheckCloseEnough(modelOut[['test44b']]$result, mxEval(omxNot(omxNot(Z)), model), 0.001)
omxCheckCloseEnough(modelOut[['test44c']]$result, mxEval(omxNot(t(Z)), model), 0.001)
omxCheckCloseEnough(modelOut[['test45']]$result, mxEval(omxSelectRows(O,N), model), 0.001)
omxCheckCloseEnough(modelOut[['test46']]$result, mxEval(omxSelectCols(O,N), model), 0.001)
omxCheckCloseEnough(modelOut[['test47']]$result, mxEval(omxSelectRowsAndCols(O,N), model), 0.001)
omxCheckCloseEnough(modelOut[['test45b']]$result, mxEval(omxSelectRows(O,omxNot(N)), model), 0.001)
omxCheckCloseEnough(modelOut[['test46b']]$result, mxEval(omxSelectCols(O,omxNot(N)), model), 0.001)
omxCheckCloseEnough(modelOut[['test47b']]$result, mxEval(omxSelectRowsAndCols(O,omxNot(N)), model), 0.001)
omxCheckCloseEnough(modelOut[['test45c']]$result, mxEval(omxSelectRows(O,N%x%5), model), 0.001)
omxCheckCloseEnough(modelOut[['test46c']]$result, mxEval(omxSelectCols(O,N%x%5), model), 0.001)
omxCheckCloseEnough(modelOut[['test47c']]$result, mxEval(omxSelectRowsAndCols(O,N%x%5), model), 0.001)
omxCheckCloseEnough(modelOut[['test45d']]$result, mxEval(omxSelectRows(O,t(N)), model), 0.001)
omxCheckCloseEnough(modelOut[['test46d']]$result, mxEval(omxSelectCols(O,t(N)), model), 0.001)
omxCheckCloseEnough(modelOut[['test47d']]$result, mxEval(omxSelectRowsAndCols(O,t(N)), model), 0.001)
omxCheckCloseEnough(modelOut[['test45e']]$result, mxEval(omxSelectRows(O,omxNot(t(N))), model), 0.001)
omxCheckCloseEnough(modelOut[['test46e']]$result, mxEval(omxSelectCols(O,omxNot(t(N))), model), 0.001)
omxCheckCloseEnough(modelOut[['test47e']]$result, mxEval(omxSelectRowsAndCols(O,omxNot(t(N))), model), 0.001)
omxCheckCloseEnough(modelOut[['test45f']]$result, mxEval(omxSelectRows(O,t(N)%x%5), model), 0.001)
omxCheckCloseEnough(modelOut[['test46f']]$result, mxEval(omxSelectCols(O,t(N)%x%5), model), 0.001)
omxCheckCloseEnough(modelOut[['test47f']]$result, mxEval(omxSelectRowsAndCols(O,t(N)%x%5), model), 0.001)
omxCheckCloseEnough(modelOut[['test45g']]$result, mxEval(omxSelectRows(O[,1:2],N%x%5), model), 0.001)
omxCheckCloseEnough(modelOut[['test46g']]$result, mxEval(omxSelectCols(O[1:2,],N%x%5), model), 0.001)
omxCheckCloseEnough(modelOut[['test45h']]$result, mxEval(omxSelectRows(O[,1:2],t(N)%x%5), model), 0.001)
omxCheckCloseEnough(modelOut[['test46h']]$result, mxEval(omxSelectCols(O[1:2,],t(N)%x%5), model), 0.001)
omxCheckCloseEnough(modelOut[['test48a']]$result, mxEval(mean(A), model), 0.001)
omxCheckCloseEnough(modelOut[['test48b']]$result, mxEval(mean(B), model), 0.001)
omxCheckCloseEnough(modelOut[['test48c']]$result, mxEval(mean(C), model), 0.001)
omxCheckCloseEnough(modelOut[['test48d']]$result, mxEval(mean(t(C)), model), 0.001)
omxCheckCloseEnough(modelOut[['test48e']]$result, mxEval(mean(-C), model), 0.001)
omxCheckCloseEnough(modelOut[['test49a']]$result, mxEval(omxGreaterThan(A,B), model), 0.001)
omxCheckCloseEnough(modelOut[['test49b']]$result, mxEval(omxGreaterThan(B,A), model), 0.001)
omxCheckCloseEnough(modelOut[['test49c']]$result, mxEval(omxGreaterThan(A,A), model), 0.001)
omxCheckCloseEnough(modelOut[['test49d']]$result, mxEval(omxGreaterThan(B,B), model), 0.001)
omxCheckCloseEnough(modelOut[['test50a']]$result, mxEval(omxLessThan(A,B), model), 0.001)
omxCheckCloseEnough(modelOut[['test50b']]$result, mxEval(omxLessThan(B,A), model), 0.001)
omxCheckCloseEnough(modelOut[['test50c']]$result, mxEval(omxLessThan(A,A), model), 0.001)
omxCheckCloseEnough(modelOut[['test50d']]$result, mxEval(omxLessThan(B,B), model), 0.001)
omxCheckCloseEnough(modelOut[['test51a']]$result, mxEval(omxAnd(omxNot(A),B), model), 0.001)
omxCheckCloseEnough(modelOut[['test51b']]$result, mxEval(omxAnd(B,A), model), 0.001)
omxCheckCloseEnough(modelOut[['test51c']]$result, mxEval(omxAnd(omxNot(A),A), model), 0.001)
omxCheckCloseEnough(modelOut[['test51d']]$result, mxEval(omxAnd(A,omxNot(A)), model), 0.001)
omxCheckCloseEnough(modelOut[['test52a']]$result, mxEval(omxOr(omxNot(A),B), model), 0.001)
omxCheckCloseEnough(modelOut[['test52b']]$result, mxEval(omxOr(B,A), model), 0.001)
omxCheckCloseEnough(modelOut[['test52c']]$result, mxEval(omxOr(omxNot(A),omxNot(A)), model), 0.001)
omxCheckCloseEnough(modelOut[['test52d']]$result, mxEval(omxOr(A,A), model), 0.001)
omxCheckCloseEnough(modelOut[['test53a']]$result, mxEval(omxApproxEquals(A,B, EPSILON), model), 0.001)
omxCheckCloseEnough(modelOut[['test53b']]$result, mxEval(omxApproxEquals(B,A, EPSILON), model), 0.001)
omxCheckCloseEnough(modelOut[['test53c']]$result, mxEval(omxApproxEquals(A,A, EPSILON), model), 0.001)
omxCheckCloseEnough(modelOut[['test53d']]$result, mxEval(omxApproxEquals(B,B, EPSILON), model), 0.001)
omxCheckCloseEnough(modelOut[['test53e']]$result, mxEval(omxApproxEquals(A,B,t(EPSILON)), model), 0.001)
omxCheckCloseEnough(modelOut[['test55a']]$result, mxEval(omxExponential(A), model), 0.001)
omxCheckCloseEnough(modelOut[['test55b']]$result, mxEval(omxExponential(t(A)), model), 0.001)
omxCheckCloseEnough(modelOut[['test55c']]$result, mxEval(omxExponential(B), model), 0.001)
omxCheckCloseEnough(modelOut[['test55d']]$result, mxEval(omxExponential(t(B)), model), 0.001)
omxCheckCloseEnough(modelOut[['test55e']]$result, mxEval(omxExponential(I), model), 0.001)
omxCheckCloseEnough(modelOut[['test55f']]$result, mxEval(omxExponential(t(I)), model), 0.001)
omxCheckCloseEnough(modelOut[['test55a']]$result, mxEval(as.matrix(expm(A)), model), 0.001)
omxCheckCloseEnough(modelOut[['test55b']]$result, mxEval(as.matrix(expm(t(A))), model), 0.001)
omxCheckCloseEnough(modelOut[['test55c']]$result, mxEval(as.matrix(expm(B)), model), 0.001)
omxCheckCloseEnough(modelOut[['test55d']]$result, mxEval(as.matrix(expm(t(B))), model), 0.001)
omxCheckCloseEnough(modelOut[['test55e']]$result, mxEval(as.matrix(expm(I)), model), 0.001)
omxCheckCloseEnough(modelOut[['test55f']]$result, mxEval(as.matrix(expm(t(I))), model), 0.001)
omxCheckCloseEnough(modelOut[['test56']]$result, mxEval(as.matrix(chol(V)), model), 0.001)
omxCheckCloseEnough(modelOut[['test57a']]$result, mxEval(as.matrix(cov2cor(A)), model), 0.001)
omxCheckCloseEnough(modelOut[['test57b']]$result, mxEval(as.matrix(cov2cor(B)), model), 0.001)
omxCheckCloseEnough(modelOut[['test58a']]$result, mxEval(vech2full(vech( t(A) %*% A )), model), 0.001)
omxCheckCloseEnough(modelOut[['test58b']]$result, mxEval(vechs2full(vechs(cov2cor(t(A) %*% A))), model), 0.001)
omxCheckCloseEnough(modelOut[['test65']]$result, mxEval(omxSelectRows(t(O),N), model), 0.001)
omxCheckCloseEnough(modelOut[['test66']]$result, mxEval(omxSelectCols(t(O),N), model), 0.001)
omxCheckCloseEnough(modelOut[['test67']]$result, mxEval(omxSelectRowsAndCols(t(O),N), model), 0.001)
omxCheckCloseEnough(modelOut[['test65b']]$result, mxEval(omxSelectRows(t(O),omxNot(N)), model), 0.001)
omxCheckCloseEnough(modelOut[['test66b']]$result, mxEval(omxSelectCols(t(O),omxNot(N)), model), 0.001)
omxCheckCloseEnough(modelOut[['test67b']]$result, mxEval(omxSelectRowsAndCols(t(O),omxNot(N)), model), 0.001)
omxCheckCloseEnough(modelOut[['test65c']]$result, mxEval(omxSelectRows(t(O),N%x%5), model), 0.001)
omxCheckCloseEnough(modelOut[['test66c']]$result, mxEval(omxSelectCols(t(O),N%x%5), model), 0.001)
omxCheckCloseEnough(modelOut[['test67c']]$result, mxEval(omxSelectRowsAndCols(t(O),N%x%5), model), 0.001)
omxCheckCloseEnough(modelOut[['test65d']]$result, mxEval(omxSelectRows(t(O),t(N)), model), 0.001)
omxCheckCloseEnough(modelOut[['test66d']]$result, mxEval(omxSelectCols(t(O),t(N)), model), 0.001)
omxCheckCloseEnough(modelOut[['test67d']]$result, mxEval(omxSelectRowsAndCols(t(O),t(N)), model), 0.001)
omxCheckCloseEnough(modelOut[['test65e']]$result, mxEval(omxSelectRows(t(O),omxNot(t(N))), model), 0.001)
omxCheckCloseEnough(modelOut[['test66e']]$result, mxEval(omxSelectCols(t(O),omxNot(t(N))), model), 0.001)
omxCheckCloseEnough(modelOut[['test67e']]$result, mxEval(omxSelectRowsAndCols(t(O),omxNot(t(N))), model), 0.001)
omxCheckCloseEnough(modelOut[['test65f']]$result, mxEval(omxSelectRows(t(O),t(N)%x%5), model), 0.001)
omxCheckCloseEnough(modelOut[['test66f']]$result, mxEval(omxSelectCols(t(O),t(N)%x%5), model), 0.001)
omxCheckCloseEnough(modelOut[['test67f']]$result, mxEval(omxSelectRowsAndCols(t(O),t(N)%x%5), model), 0.001)
omxCheckCloseEnough(modelOut[['test68a']]$result, mxEval(logm(A), model), .001)
omxCheckCloseEnough(modelOut[['test68b']]$result, mxEval(logm(B), model), .001)
omxCheckCloseEnough(modelOut[['test68c']]$result, mxEval(logm(I), model), .001)
omxCheckCloseEnough(modelOut[['test69a']]$result, asin(A$values), .001)
omxCheckCloseEnough(modelOut[['test69b']]$result, acos(A$values), .001)
omxCheckCloseEnough(modelOut[['test69c']]$result, atan(A$values), .001)
omxCheckCloseEnough(modelOut[['test70']]$result, lgamma(A$values), .001)
omxCheckCloseEnough(modelOut[['test71a']]$result, qnorm(p=(A$values)), .001)
omxCheckCloseEnough(modelOut[['test72a']]$result, asinh(A$values), .001)
omxCheckCloseEnough(modelOut[['test72b']]$result, acosh(D$values), .001)
omxCheckCloseEnough(modelOut[['test72c']]$result, atanh(A$values), .001)
omxCheckCloseEnough(modelOut[['test73']]$result, lgamma(A$values+1), .001)
omxCheckCloseEnough(modelOut[['test71b']]$result, qnorm(log(A$values),log.p=T), .001)
omxCheckCloseEnough(modelOut[['test71b']]$result, mxEval(logp2z(log(A)),model,T), .001)
omxCheckCloseEnough(modelOut[['test74a']]$result, dbeta(x=A$values,shape1=A$values,shape2=A$values,ncp=A$values,log=T),
										.001)
omxCheckCloseEnough(modelOut[['test74a']]$result, mxEval(dbeta(A,A,A,A,T),model,T), .001)
omxCheckCloseEnough(modelOut[['test75a']]$result, pbeta(q=A$values,shape1=A$values,shape2=A$values,ncp=A$values,
																										lower.tail=T,log.p=T),.001)
omxCheckCloseEnough(modelOut[['test75a']]$result, mxEval(pbeta(A,A,A,A,T,T),model,T), .001)
omxCheckCloseEnough(modelOut[['test74b']]$result, dbeta(A$values,1,1), .001)
omxCheckCloseEnough(modelOut[['test75b']]$result, pbeta(A$values,1,1), .001)
omxCheckCloseEnough(modelOut[['test74c']]$result, dbeta(C$values,D$values,D$values), .001)
omxCheckCloseEnough(modelOut[['test75c']]$result, pbeta(C$values,D$values,D$values), .001)
omxCheckCloseEnough(modelOut[['test76a']]$result, besselI(A$values,A$values,A$values), .001)
omxCheckCloseEnough(modelOut[['test76b']]$result, besselI(C$values,D$values,0), .001)
omxCheckCloseEnough(modelOut[['test77a']]$result, besselJ(A$values,A$values), .001)
omxCheckCloseEnough(modelOut[['test77b']]$result, besselJ(C$values,D$values), .001)
omxCheckCloseEnough(modelOut[['test78a']]$result, besselK(A$values,A$values,A$values), .001)
omxCheckCloseEnough(modelOut[['test78b']]$result, besselK(C$values,D$values,0), .001)
omxCheckCloseEnough(modelOut[['test79a']]$result, besselY(A$values,A$values), .001)
omxCheckCloseEnough(modelOut[['test79b']]$result, besselY(C$values,D$values), .004)
omxCheckCloseEnough(modelOut[['test80a']]$result, dpois(D$values,D$values), .001)
omxCheckCloseEnough(modelOut[['test80b']]$result, dpois(D$values,D$values,TRUE), .001)
omxCheckCloseEnough(modelOut[['test81a']]$result, ppois(D$values,D$values,TRUE,TRUE), .001)
omxCheckCloseEnough(modelOut[['test81b']]$result, ppois(D$values,D$values,TRUE,FALSE), .001)
omxCheckCloseEnough(modelOut[['test81c']]$result, ppois(D$values,D$values,FALSE,FALSE), .001)
omxCheckCloseEnough(modelOut[['test82a']]$result, dnbinom(x=D$values,size=D$values,prob=C$values[1:10],log=FALSE), .001)
omxCheckCloseEnough(modelOut[['test82b']]$result, dnbinom(x=D$values,size=D$values,prob=C$values[1:10],log=TRUE), .001)
omxCheckCloseEnough(modelOut[['test83a']]$result, pnbinom(q=D$values,size=D$values,prob=C$values[1:10],lower.tail=TRUE,
																											 log.p=FALSE), .001)
omxCheckCloseEnough(modelOut[['test83b']]$result, pnbinom(q=D$values,size=D$values,prob=C$values[1:10],lower.tail=TRUE,
																											 log.p=FALSE), .001)
omxCheckCloseEnough(modelOut[['test83c']]$result, pnbinom(q=D$values,size=D$values,prob=C$values[1:10],lower.tail=FALSE,
																											 log.p=TRUE), .001)
omxCheckCloseEnough(modelOut[['test84a']]$result, dnbinom(x=D$values,size=D$values,mu=C$values[1:10],log=FALSE), .001)
omxCheckCloseEnough(modelOut[['test84b']]$result, dnbinom(x=D$values,size=D$values,mu=C$values[1:10],log=TRUE), .001)
omxCheckCloseEnough(modelOut[['test85a']]$result, pnbinom(q=D$values,size=D$values,mu=C$values[1:10],lower.tail=TRUE,
																											 log.p=FALSE), .001)
omxCheckCloseEnough(modelOut[['test85b']]$result, pnbinom(q=D$values,size=D$values,mu=C$values[1:10],lower.tail=TRUE,
																											 log.p=FALSE), .001)
omxCheckCloseEnough(modelOut[['test85c']]$result, pnbinom(q=D$values,size=D$values,mu=C$values[1:10],lower.tail=FALSE,
																											 log.p=TRUE), .001)
omxCheckCloseEnough(modelOut[['test86a']]$result, dchisq(x=A$values,df=A$values,ncp=A$values,log=TRUE), .001)
omxCheckCloseEnough(modelOut[['test86b']]$result, pchisq(q=A$values,df=A$values,ncp=A$values,lower.tail=TRUE,log.p=TRUE), .001)
omxCheckCloseEnough(modelOut[['test87a']]$result, dchisq(x=D$values,df=1), .001)
omxCheckCloseEnough(modelOut[['test87b']]$result, pchisq(q=D$values,df=1,lower.tail=FALSE), .001)
omxCheckCloseEnough(modelOut[['test88a']]$result, dbinom(x=D$values,size=D$values,prob=A$values[1:10],log=FALSE), .001)
omxCheckCloseEnough(modelOut[['test88b']]$result, dbinom(x=D$values,size=D$values,prob=A$values[1:10],log=TRUE), .001)
omxCheckCloseEnough(modelOut[['test89a']]$result, pbinom(q=D$values,size=D$values,prob=A$values[1:10],lower.tail=FALSE,log.p=FALSE), .001)
omxCheckCloseEnough(modelOut[['test89b']]$result, pbinom(q=D$values,size=D$values,prob=A$values[1:10],lower.tail=TRUE,log.p=TRUE), .001)
omxCheckCloseEnough(modelOut[['test90a']]$result, dcauchy(x=H$values,location=A$values[1:16],scale=D$values[c(1:10,1:6)],log=FALSE), .001)
omxCheckCloseEnough(modelOut[['test90b']]$result, dcauchy(x=H$values,location=A$values[1:16],scale=D$values[c(1:10,1:6)],log=TRUE), .001)
omxCheckCloseEnough(modelOut[['test91a']]$result, pcauchy(q=H$values,location=A$values[1:16],scale=D$values[c(1:10,1:6)],lower.tail=FALSE,log=TRUE),
										.001)
omxCheckCloseEnough(modelOut[['test91b']]$result, pcauchy(q=H$values,location=A$values[1:16],scale=D$values[c(1:10,1:6)],lower.tail=TRUE,log=FALSE),
										.001)


#Test for regressions in compatibility between mxEval() and negative-binomial probability functions:
omxCheckCloseEnough(modelOut[["test82a"]]$result, mxEval(omxDnbinom(D,D,C,-1,0),modelOut,T))
omxCheckCloseEnough(modelOut[["test82b"]]$result, mxEval(omxDnbinom(D,D,C,-1,C),modelOut,T))
omxCheckCloseEnough(modelOut[["test83a"]]$result, mxEval(omxPnbinom(D,D,C,-1,C,0),modelOut,T))
omxCheckCloseEnough(modelOut[["test83b"]]$result, mxEval(omxPnbinom(D,D,C,-1,1,0),modelOut,T))
omxCheckCloseEnough(modelOut[["test83c"]]$result, mxEval(omxPnbinom(D,D,C,-1,0,1),modelOut,T))
omxCheckCloseEnough(modelOut[["test84a"]]$result, mxEval(omxDnbinom(D,D,-1,C,0),modelOut,T))
omxCheckCloseEnough(modelOut[["test84b"]]$result, mxEval(omxDnbinom(D,D,-1,C,C),modelOut,T))
omxCheckCloseEnough(modelOut[["test85a"]]$result, mxEval(omxPnbinom(D,D,-1,C,C,0),modelOut,T))
omxCheckCloseEnough(modelOut[["test85b"]]$result, mxEval(omxPnbinom(D,D,-1,C,1,0),modelOut,T))
omxCheckCloseEnough(modelOut[["test85c"]]$result, mxEval(omxPnbinom(D,D,-1,C,0,1),modelOut,T))
omxCheckCloseEnough(modelOut[["test92"]]$result, mxEval(rowSums(C),modelOut,T))
omxCheckCloseEnough(modelOut[["test93"]]$result, mxEval(colSums(C),modelOut,T))

# Test for allowable scalar multiplication in frontend and backend.
omxCheckCloseEnough(modelOut[['test94']]$result, mxEval(test94,model,T), .001)
omxCheckCloseEnough(modelOut[['test94']]$result, 1/(1+N$values), 0.001)


# Check internal function for definition variables
m0 <- mxModel()
omxCheckTrue(!imxHasDefinitionVariable(m0))
m1 <- mxModel(mxMatrix('Full', 1, 1, labels='data.blah'))
omxCheckTrue(imxHasDefinitionVariable(m1)) # still has no data
m2 <- mxModel(m1, mxData(matrix(0, 1, 1), 'raw'))
omxCheckTrue(imxHasDefinitionVariable(m2))
m3 <- mxModel(name='b', m1)
omxCheckTrue(imxHasDefinitionVariable(m3)) # submodel behavior
m3 <- mxModel(name='b', m2)
omxCheckTrue(imxHasDefinitionVariable(m3))
m4 <- mxModel(mxAlgebra(data.x1 %x% V, name='bob'), mxData(matrix(0, 1, 1), 'raw'))
omxCheckTrue(imxHasDefinitionVariable(m4)) # algebra

# ensure model$output is correctly populated
omxCheckCloseEnough(dim(modelOut$output$algebras[['model.test0']]), c(5,5), 1)