Added license information to R source files.
[openmx:openmx.git] / R / MxMatrix.R
1 #
2 #   Copyright 2007-2009 The OpenMx Project
3 #
4 #   Licensed under the Apache License, Version 2.0 (the "License");
5 #   you may not use this file except in compliance with the License.
6 #   You may obtain a copy of the License at
7
8 #        http://www.apache.org/licenses/LICENSE-2.0
9
10 #    Unless required by applicable law or agreed to in writing, software
11 #    distributed under the License is distributed on an "AS IS" BASIS,
12 #    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 #   See the License for the specific language governing permissions and
14 #   limitations under the License.
15
16
17 single.na <- function(a) {\r
18         return((length(a) == 1) && \r
19                 (is.list(a) || is.vector(a)) && \r
20                 (is.na(a) == TRUE))\r
21 }\r
22 \r
23 is.Matrix <- function(a) {\r
24         return(is.matrix(a) || is(a, "Matrix")) \r
25 }\r
26 \r
27 verifySquare <- function(.Object) {\r
28         if (nrow(.Object@specification) != ncol(.Object@specification)) { 
29                 stop(paste("Specification matrix of MxMatrix", 
30                                 omxQuotes(.Object@name), "is not square"), call.=FALSE)
31         }\r
32         if (nrow(.Object@values) != ncol(.Object@values)) {
33                 stop(paste("Values matrix of MxMatrix", 
34                                 omxQuotes(.Object@name), "is not square"), call.=FALSE)
35         }               \r
36 }\r
37 \r
38 setGeneric("omxVerifyMatrix", function(.Object) { \r
39         return(standardGeneric("omxVerifyMatrix")) \r
40 } )\r
41 \r
42 #\r
43 # The Matrix package returns a non-symmetric matrix\r
44 # when you modify a symmetric matrix.\r
45 # This prevents us from using symmetric matrices from the package.\r
46 #\r
47 setClass(Class = "MxSymmetricMatrix",\r
48         representation = representation(\r
49                 specification = "MxSymmetricSparse",\r
50                 values = "MxSymmetricSparse", name = "character", "VIRTUAL"))\r
51 \r
52 setClass(Class = "MxNonSymmetricMatrix",\r
53         representation = representation(\r
54                 specification = "MxSparseMatrix",\r
55                 values = "Matrix", name = "character", "VIRTUAL"))\r
56                 \r
57 setClassUnion("MxMatrix",\r
58     c("MxSymmetricMatrix", "MxNonSymmetricMatrix"))\r
59                 \r
60                 \r
61 setMethod("initialize", "MxMatrix",\r
62         function(.Object, specification, values, name) {\r
63                 .Object@specification = specification\r
64                 .Object@values = values\r
65                 .Object@name = name\r
66                 return(.Object)\r
67         }\r
68 )               \r
69 \r
70 setMethod("omxVerifyMatrix", "MxMatrix",\r
71         function(.Object) {\r
72                 if (nrow(.Object@specification) != nrow(.Object@values)) {\r
73                         stop(paste("Specification and values matrices of", \r
74                                 omxQuotes(.Object@name), "have different dimensions"), call.=FALSE)\r
75                 }\r
76                 if (ncol(.Object@specification) != ncol(.Object@values)) {\r
77                         stop(paste("Specification and values matrices of", \r
78                                 omxQuotes(.Object@name), "have different dimensions"), call.=FALSE)\r
79                 }               \r
80         }\r
81 )\r
82 \r
83 setMethod("nrow", "MxMatrix",\r
84         function(x) {\r
85             return(nrow(x@specification))\r
86         }\r
87 )\r
88 \r
89 setMethod("ncol", "MxMatrix",\r
90         function(x) {\r
91             return(ncol(x@specification))\r
92         }\r
93 )\r
94 \r
95 setMethod("[", "MxMatrix",\r
96         function(x, i, j, ..., drop = FALSE) {\r
97                 return(x@values[i,j])\r
98     }\r
99 )\r
100 \r
101 setReplaceMethod("[", "MxMatrix", \r
102         function(x, i, j, value) {\r
103                 x@values[i,j] <- value\r
104                 return(x)\r
105     }\r
106 )\r
107 \r
108 matrixTypes <- c("Diag", "Full", "Iden", "Symm", "Unit", "Zero")\r
109 squareMatrices <- c("Diag", "Iden", "Symm")\r
110 \r
111 \r
112 mxMatrix <- function(type = "Full", values = NA, \r
113         specification = NA, nrow = NA, ncol = NA,\r
114         byrow = FALSE, free = FALSE, name = NA) {\r
115         if (is.na(match(type, matrixTypes))) {\r
116                 stop(paste("Type must be one of:", paste(matrixTypes, collapse=" ")))\r
117         }\r
118         if ((is.Matrix(values) || is.Matrix(specification)) &&\r
119                 (!is.null(match.call()$nrow))) {\r
120                 warning("\'nrow\' is disregarded for mxMatrix constructor")\r
121         }\r
122         if ((is.Matrix(values) || is.Matrix(specification)) &&\r
123                 (!is.null(match.call()$ncol))) {\r
124                 warning("\'ncol\' is disregarded for mxMatrix constructor")\r
125         }\r
126         if (is.Matrix(values) && is.Matrix(specification) &&\r
127                 !all(dim(values) == dim(specification))) {\r
128                 stop("Values and specification matrices are not of identical dimensions")\r
129         }\r
130         if (is.Matrix(values)) {\r
131                 nrow <- nrow(values)\r
132                 ncol <- ncol(values)\r
133         } else if (is.Matrix(specification)) {\r
134                 nrow <- nrow(specification)\r
135                 ncol <- ncol(specification)\r
136         }\r
137         if (type %in% squareMatrices) {\r
138                 if (is.na(nrow) && is.na(ncol)) {\r
139                         stop("Either nrow or ncol must be specified on a square matrix")\r
140                 } else if (is.na(nrow)) {\r
141                         nrow <- ncol\r
142                 } else if (is.na(ncol)) {\r
143                         ncol <- nrow\r
144                 }\r
145         } else if (is.na(nrow) || is.na(ncol)) {\r
146                 stop("Both nrow and ncol must be specified on a non-square matrix")\r
147         }
148         if (byrow) {
149                 values <- omxCreateByRow(values, nrow, ncol, type)
150                 specification <- omxCreateByRow(specification, nrow, ncol, type)\r
151         }\r
152         if (is.na(name)) {\r
153                 name <- omxUntitledName()\r
154         }\r
155         if (!is.character(name)) {\r
156                 stop("\'name\' must be a character vector!")\r
157         }\r
158         typeName <- paste(type, "Matrix", sep="")\r
159         return(new(typeName, name, values, specification, \r
160                         nrow, ncol, byrow, free))\r
161 }\r
162
163 omxCreateByRow <- function(source, nrow, ncol, type) {
164         if (!is.na(source) && is.vector(source)) {
165                 if (type %in% squareMatrices) {
166                         res <- matrix(0, nrow, ncol)
167                         res[upper.tri(res, diag=TRUE)] <- source
168                         res[lower.tri(res)] <- res[upper.tri(res)]
169                         return(res)
170                 } else {
171                         return(matrix(source, nrow, ncol, TRUE))
172                 }
173         } else {
174                 return(source)
175         }
176 }\r
177 \r
178 omxSparseMatrixParameters <- function(specification, bounds, result, 
179         defNames, matrixNumber, reverse=FALSE) {
180         if (length(specification@dataVector) == 0) {\r
181                 return(result)\r
182         }\r
183         for(i in 1:length(specification@dataVector)) {\r
184             if (reverse == FALSE || specification@rowVector[i] != 
185                 specification@colVector[i]) {\r
186                     parameterName <- as.character(specification@dataVector[[i]])\r
187                     if (reverse) {\r
188                             col <- specification@rowVector[i] - 1\r
189                             row <- specification@colVector[i] - 1\r
190                         } else {\r
191                             row <- specification@rowVector[i] - 1\r
192                             col <- specification@colVector[i] - 1                       \r
193                         }\r
194                         boundsLookup <- omxLocateBounds(bounds, parameterName)\r
195                         minBounds <- boundsLookup[[1]]\r
196                         maxBounds <- boundsLookup[[2]]\r
197                         if (is.na(parameterName)) {\r
198                                 result[[length(result)+1]] <-  list(minBounds, maxBounds, \r
199                                         c(matrixNumber, row, col))\r
200                         } else if (!(parameterName %in% defNames)) {\r
201                                 if (!is.null(result[[parameterName]])) {\r
202                                         original <- result[[parameterName]]\r
203                                         original[[length(original) + 1]] <- c(matrixNumber, row, col)\r
204                                         result[[parameterName]] <- original\r
205                                 } else {\r
206                                         result[[parameterName]] <- list(minBounds, maxBounds, \r
207                                                 c(matrixNumber, row,col))\r
208                                 }\r
209                         }\r
210                 }\r
211         }\r
212         return(result)\r
213 }\r
214
215 omxSparseMatrixDefinitions <- function(specification, bounds, result, 
216         defLocations, matrixNumber, reverse=FALSE) {
217         if (length(specification@dataVector) == 0) {
218                 return(result)
219         }
220         defNames <- names(defLocations)
221         for(i in 1:length(specification@dataVector)) {
222             if (reverse == FALSE || specification@rowVector[i] != 
223                 specification@colVector[i]) {
224                     parameterName <- as.character(specification@dataVector[[i]])
225                     if (reverse) {
226                             col <- specification@rowVector[i] - 1
227                             row <- specification@colVector[i] - 1
228                         } else {
229                             row <- specification@rowVector[i] - 1
230                             col <- specification@colVector[i] - 1                       
231                         }
232                         if (parameterName %in% defNames) {
233                                 if (!is.null(result[[parameterName]])) {
234                                         original <- result[[parameterName]]
235                                         dataNumber <- 
236                                         original[[length(original) + 1]] <- c(matrixNumber, row, col)
237                                         result[[parameterName]] <- original
238                                 } else {
239                                         dataNumber <- defLocations[[parameterName]][[1]]
240                                         columnNumber <- defLocations[[parameterName]][[2]]
241                                         result[[parameterName]] <- list(dataNumber, columnNumber, 
242                                                 c(matrixNumber, row,col))
243                                 }
244                         }
245                 }
246         }
247         return(result)
248 }
249 \r
250 generateMatrixListHelper <- function(mxMatrix) {\r
251         values <- mxMatrix@values\r
252     if(is(values,"MxSymmetricSparse")) {\r
253                 return(as.matrix(values))    #return(Matrix(as.matrix(values)))    \r
254     } else {\r
255                 return(as.matrix(values))\r
256     }\r
257 }\r
258 \r
259 omxGenerateParameterListHelper <- function(mxMatrix, bounds,\r
260         result, defNames, matrixNumber) {\r
261         specification <- mxMatrix@specification\r
262         result <- omxSparseMatrixParameters(specification, bounds,\r
263                 result, defNames, matrixNumber)\r
264         if(is(specification,"MxSymmetricSparse")) {\r
265                 result <- omxSparseMatrixParameters(specification, bounds, \r
266                         result, defNames, matrixNumber, TRUE)\r
267         }\r
268         return(result)\r
269 }\r
270
271 omxGenerateDefinitionListHelper <- function(mxMatrix, bounds,
272         result, defLocations, matrixNumber) {
273         specification <- mxMatrix@specification
274         result <- omxSparseMatrixDefinitions(specification, bounds,
275                 result, defLocations, matrixNumber)
276         if(is(specification,"MxSymmetricSparse")) {
277                 result <- omxSparseMatrixDefinitions(specification, bounds, 
278                         result, defLocations, matrixNumber, TRUE)
279         }
280         return(result)
281 }
282 \r
283 omxDisplayMatrix <- function(mxMatrix) {\r
284    cat("MxMatrix", omxQuotes(mxMatrix@name), '\n')
285    cat("\n")\r
286    cat("Specification matrix:\n")\r
287    print(mxMatrix@specification, use.quotes = TRUE)
288    cat("\n")\r
289    cat("Values matrix:\n")\r
290    values <- mxMatrix@values\r
291    if(is(values, "sparseMatrix")) {\r
292       print(as(values, 'matrix'))\r
293    } else if(is(values, "MxSymmetricSparse")) {\r
294       print(values, use.quotes = FALSE)\r
295    } else {\r
296       print(values)\r
297    }\r
298    cat("\n")\r
299 }\r
300 \r
301 setMethod("print", "MxMatrix", function(x,...) { omxDisplayMatrix(x) })\r
302 setMethod("show", "MxMatrix", function(object) { omxDisplayMatrix(object) })\r