1### SUMT (Sequential Unconstrained Maximization Technique)
2### borrowed from package 'clue'
3###
4### Adapted for linear constraints
5sumt <- function(fn, grad=NULL, hess=NULL,
6                 start,
7                 maxRoutine, constraints,
8                 SUMTTol = sqrt(.Machine$double.eps),
9                           # difference between estimates for successive outer iterations
10                 SUMTPenaltyTol = sqrt(.Machine$double.eps),
11                           # maximum allowed penalty
12                 SUMTQ = 10,
13                 SUMTRho0 = NULL,
14                 printLevel=print.level, print.level=0,
15                 SUMTMaxIter=100,
16                 ...) {
17   ## constraints    list w/components eqA and eqB.  Maximization will
18   ##                be performed wrt to the constraint
19   ##                A %*% theta + B = 0
20   ##                The user must ensure the matrices are in correct
21   ##                form
22   ## maxSUMTiter    how many SUMT iterations to perform max
23   ##
24   penalty <- function(theta) {
25      p <- A %*% theta + B
26      sum(p*p)
27   }
28   ## Penalty gradient and Hessian are used only if corresponding function
29   ## for the likelihood function is provided
30   gPenalty <- function(theta) {
31      2*(t(theta) %*% t(A) %*% A - t(B) %*% A)
32   }
33   hessPenalty <- function(theta) {
34      2*t(A) %*% A
35   }
36
37   ## strip possible arguments of maxRoutine and call the function thereafter
38   callWithoutMaxArgs <- function(theta, fName, ...) {
39      return( callWithoutArgs( theta, fName = fName,
40         args = names(formals(maxRoutine)), ... ) )
41   }
42
43   SUMTMessage <- function(code) {
44      message <- switch(code,
45                        "1" = "penalty close to zero",
46                        "2" = "successive function values within tolerance limit",
47                        "4" = "Outer iteration limit exceeded (increase SUMTMaxIter ?).",
48                        paste("Code", code))
49      return(message)
50   }
51   ## the penalized objective function
52   Phi <- function(theta, ...) {
53      llVal <- callWithoutMaxArgs( theta, "logLikFunc", fnOrig = fn,
54         gradOrig = grad, hessOrig = hess, sumObs = FALSE, ... )
55      llVal <- llVal - rho * penalty( theta ) / length( llVal )
56      g <- attributes( llVal )$gradient
57      if( !is.null( g ) ) {
58         if( is.matrix( g ) ) {
59            g <- g - matrix(
60               rep( rho * gPenalty( theta ) / nrow( g ), each = nrow( g ) ),
61               nrow = nrow( g ), ncol = ncol( g ) )
62         } else {
63            g <- g - rho * gPenalty( theta )
64         }
65         attributes( llVal )$gradient <- g
66      }
67      h <- attributes( llVal )$hessian
68      if( !is.null( h ) ) {
69         attributes( llVal )$hessian <- h - rho * hessPenalty( theta )
70      }
71      return( llVal )
72   }
73   ## gradient of the penalized objective function
74   if(!is.null(grad)) {
75      gradPhi<- function(theta, ...) {
76         g <- grad(theta, ...)
77         if(is.matrix(g)) {
78            g <- g - matrix(
79               rep( rho * gPenalty( theta ) / nrow( g ), each = nrow( g ) ),
80               nrow = nrow( g ), ncol = ncol( g ) )
81         } else {
82            g <- g - rho * gPenalty( theta )
83         }
84         return( g )
85      }
86   } else {
87      gradPhi <- NULL
88   }
89   ## Hessian of the penalized objective function
90   if(!is.null(hess)) {
91      hessPhi <- function(theta, ...) {
92         return( hess(theta, ...) - rho*hessPenalty(theta) )
93      }
94   } else {
95      hessPhi <- NULL
96   }
97   ## -------- SUMT Main code ---------
98   ## Note also that currently we do not check whether optimization was
99   ## "successful" ...
100   A <- constraints$eqA
101   B <- as.matrix(constraints$eqB)
102   ## Check if the matrices conform
103   if(ncol(A) != length(start)) {
104      stop("Equality constraint matrix A must have the same number\n",
105           "of columns as the parameter length ",
106           "(currently ", ncol(A), " and ", length(start), ")")
107   }
108   if(nrow(A) != nrow(B)) {
109      stop("Equality constraint matrix A must have the same number\n",
110           "of rows as the matrix B ",
111           "(currently ", nrow(A), " and ", nrow(B), ")")
112   }
113   ## Find a suitable inital value for rho if not specified
114   if(is.null(SUMTRho0)) {
115      rho <- 0
116      result <- maxRoutine(fn=Phi, grad=gradPhi, hess=hessPhi,
117                           start=start,
118                           printLevel=max(printLevel - 1, 0),
119                           ...)
120      theta <- coef(result)
121                           # Note: this may be a bad idea,
122                           # if unconstrained function is unbounded
123                           # from above.  In that case rather specify SUMTRho0.
124      if(printLevel > 0) {
125         cat("SUMT initial: rho = ", rho,
126             ", function = ",
127             callWithoutMaxArgs( theta, "logLikFunc",
128                                fnOrig = fn, gradOrig = grad,
129                                hessOrig = hess, ... ),
130             ", penalty = ", penalty(theta), "\n")
131         cat("Estimate:")
132         print(theta)
133      }
134      ## Better upper/lower bounds for rho?
135      rho <- max( callWithoutMaxArgs( theta, "logLikFunc", fnOrig = fn,
136         gradOrig = grad, hessOrig = hess, ... ), 1e-3) /
137         max(penalty(start), 1e-3)
138   }
139   ## if rho specified, simply pick that and use previous initial values
140   else {
141       rho <- SUMTRho0
142       theta <- start
143   }
144   ## </TODO>
145   iter <- 1L
146   repeat {
147      thetaOld <- theta
148      result <- maxRoutine(fn=Phi, grad=gradPhi, hess=hessPhi,
149                      start=thetaOld,
150                      printLevel=max(printLevel - 1, 0),
151                      ...)
152      theta <- coef(result)
153      if(printLevel > 0) {
154         cat("SUMT iteration ", iter,
155             ": rho = ", rho, ", function = ", callWithoutMaxArgs( theta,
156             "logLikFunc", fnOrig = fn, gradOrig = grad, hessOrig = hess, ... ),
157             ", penalty = ", penalty(theta), "\n", sep="")
158         cat("Estimate:")
159         print(theta)
160      }
161      if(max(abs(thetaOld - theta)) < SUMTTol) {
162         SUMTCode <- 2
163         break
164      }
165      if(penalty(theta) < SUMTPenaltyTol) {
166         SUMTCode <- 1
167         break
168      }
169      if(iter >= SUMTMaxIter) {
170         SUMTCode <- 4
171         break
172      }
173      iter <- iter + 1L
174      rho <- SUMTQ * rho
175   }
176   ## Now we replace the resulting gradient and Hessian with those,
177   ## calculated on the original function
178   llVal <- callWithoutMaxArgs( theta, "logLikFunc", fnOrig = fn,
179      gradOrig = grad, hessOrig = hess, sumObs = FALSE, ... )
180   gradient <- attr( llVal, "gradient" )
181   if( is.null( gradient ) ) {
182      gradient <- callWithoutMaxArgs( theta, "logLikGrad", fnOrig = fn,
183         gradOrig = grad, hessOrig = hess, sumObs = FALSE, ... )
184   }
185   if( !is.null( dim( gradient ) ) ) {
186      if( nrow( gradient ) > 1 ) {
187         gradientObs <- gradient
188      }
189      gradient <- colSums( gradient )
190   } else if( length( start ) == 1 && length( gradient ) > 1 ) {
191      gradientObs <- matrix( gradient, ncol = 1 )
192      gradient <- sum( gradient )
193   }
194   result$gradient <- gradient
195   names( result$gradient ) <- names( result$estimate )
196
197   result$hessian <- callWithoutMaxArgs( theta, "logLikHess", fnOrig = fn,
198      gradOrig = grad, hessOrig = hess, ... )
199   result$constraints <- list(type="SUMT",
200                             barrier.value=penalty(theta),
201                              code=SUMTCode,
202                              message=SUMTMessage(SUMTCode),
203                             outer.iterations=iter
204                             )
205   if( exists( "gradientObs" ) ) {
206      result$gradientObs <- gradientObs
207      colnames( result$gradientObs ) <- names( result$estimate )
208   }
209
210   if( result$constraints$barrier.value > 0.001 ) {
211      warning( "problem in imposing equality constraints: the constraints",
212         " are not satisfied (barrier value = ",
213         result$constraints$barrier.value, "). Try setting 'SUMTTol' to 0" )
214   }
215   return(result)
216}
217