util-SSR: SemiSmooth Reformulation
In GNE: Computation of Generalized Nash Equilibria
SSR R Documentation
SemiSmooth Reformulation
Description
functions of the SemiSmooth Reformulation of the GNEP
Usage
funSSR(z, dimx, dimlam, grobj, arggrobj, constr, argconstr, grconstr, arggrconstr,
compl, argcompl, dimmu, joint, argjoint, grjoint, arggrjoint, echo=FALSE)
jacSSR(z, dimx, dimlam, heobj, argheobj, constr, argconstr, grconstr, arggrconstr,
heconstr, argheconstr, gcompla, gcomplb, argcompl, dimmu, joint, argjoint,
grjoint, arggrjoint, hejoint, arghejoint, echo=FALSE)
Arguments
z
a numeric vector z containing (x, lambda, mu) values.
dimx
a vector of dimension for x.
dimlam
a vector of dimension for lambda.
grobj
gradient of the objective function, see details.
arggrobj
a list of additional arguments of the objective gradient.
constr
constraint function, see details.
argconstr
a list of additional arguments of the constraint function.
grconstr
gradient of the constraint function, see details.
arggrconstr
a list of additional arguments of the constraint gradient.
compl
the complementarity function with (at least) two arguments: compl(a,b).
argcompl
list of possible additional arguments for compl.
dimmu
a vector of dimension for mu.
joint
joint function, see details.
argjoint
a list of additional arguments of the joint function.
grjoint
gradient of the joint function, see details.
arggrjoint
a list of additional arguments of the joint gradient.
heobj
Hessian of the objective function, see details.
argheobj
a list of additional arguments of the objective Hessian.
heconstr
Hessian of the constraint function, see details.
argheconstr
a list of additional arguments of the constraint Hessian.
gcompla
derivative of the complementarity function w.r.t. the first argument.
gcomplb
derivative of the complementarity function w.r.t. the second argument.
hejoint
Hessian of the joint function, see details.
arghejoint
a list of additional arguments of the joint Hessian.
echo
a logical to show some traces.
Details
Compute the SemiSmooth Reformulation of the GNEP: the Generalized Nash equilibrium problem is defined
by objective functions Obj with player variables x defined in dimx and
may have player-dependent constraint functions g of dimension dimlam
and/or a common shared joint function h of dimension dimmu,
where the Lagrange multiplier are lambda and mu, respectively,
see F. Facchinei et al.(2009) where there is no joint function.
- Arguments of the Phi function
-
The arguments which are functions must respect the following features
grobj-
The gradient Grad Obj of an objective function Obj (to be minimized) must have 3 arguments for Grad Obj(z, playnum, ideriv): vector z, player number, derivative index
, and optionnally additional arguments in arggrobj.
constr-
The constraint function g must have 2 arguments: vector z, player number,
such that g(z, playnum) <= 0. Optionnally, g may have additional arguments in argconstr.
grconstr-
The gradient of the constraint function g must have 3 arguments: vector z, player number, derivative index,
and optionnally additional arguments in arggrconstr.
complIt must have two arguments and optionnally additional arguments in argcompl.
A typical example is the minimum function.
joint-
The constraint function h must have 1 argument: vector z,
such that h(z) <= 0. Optionnally, h may have additional arguments in argjoint.
grjoint-
The gradient of the constraint function h must have 2 arguments: vector z, derivative index,
and optionnally additional arguments in arggrjoint.
- Arguments of the Jacobian of Phi
-
The arguments which are functions must respect the following features
heobjIt must have 4 arguments: vector z, player number, two derivative indexes and optionnally additional arguments in argheobj.
heconstrIt must have 4 arguments: vector z, player number, two derivative indexes and optionnally additional arguments in argheconstr.
gcompla,gcomplbIt must have two arguments and optionnally additional arguments in argcompl.
hejointIt must have 3 arguments: vector z, two derivative indexes and optionnally additional arguments in arghejoint.
See the example below.
Value
A vector for funSSR or a matrix for jacSSR.
Author(s)
Christophe Dutang
References
F. Facchinei, A. Fischer and V. Piccialli (2009),
Generalized Nash equilibrium problems and Newton methods,
Math. Program.
See Also
See also GNE.nseq.
Examples
# (1) associated objective functions
#
dimx <- c(2, 2, 3)
#Gr_x_j O_i(x)
grfullob <- function(x, i, j)
{
x <- x[1:7]
if(i == 1)
{
grad <- 3*(x - 1:7)^2
}
if(i == 2)
{
grad <- 1:7*(x - 1:7)^(0:6)
}
if(i == 3)
{
s <- x[5]^2 + x[6]^2 + x[7]^2 - 5
grad <- c(1, 0, 1, 0, 4*x[5]*s, 4*x[6]*s, 4*x[7]*s)
}
grad[j]
}
#Gr_x_k Gr_x_j O_i(x)
hefullob <- function(x, i, j, k)
{
x <- x[1:7]
if(i == 1)
{
he <- diag( 6*(x - 1:7) )
}
if(i == 2)
{
he <- diag( c(0, 2, 6, 12, 20, 30, 42)*(x - 1:7)^c(0, 0:5) )
}
if(i == 3)
{
s <- x[5]^2 + x[6]^2 + x[7]^2
he <- rbind(rep(0, 7), rep(0, 7), rep(0, 7), rep(0, 7),
c(0, 0, 0, 0, 4*s+8*x[5]^2, 8*x[5]*x[6], 8*x[5]*x[7]),
c(0, 0, 0, 0, 8*x[5]*x[6], 4*s+8*x[6]^2, 8*x[6]*x[7]),
c(0, 0, 0, 0, 8*x[5]*x[7], 8*x[6]*x[7], 4*s+8*x[7]^2) )
}
he[j,k]
}
# (2) constraint linked functions
#
dimlam <- c(1, 2, 2)
#constraint function g_i(x)
g <- function(x, i)
{
x <- x[1:7]
if(i == 1)
res <- sum( x^(1:7) ) -7
if(i == 2)
res <- c(sum(x) + prod(x) - 14, 20 - sum(x))
if(i == 3)
res <- c(sum(x^2) + 1, 100 - sum(x))
res
}
#Gr_x_j g_i(x)
grfullg <- function(x, i, j)
{
x <- x[1:7]
if(i == 1)
{
grad <- (1:7) * x ^ (0:6)
}
if(i == 2)
{
grad <- 1 + sapply(1:7, function(i) prod(x[-i]))
grad <- cbind(grad, -1)
}
if(i == 3)
{
grad <- cbind(2*x, -1)
}
if(i == 1)
res <- grad[j]
if(i != 1)
res <- grad[j,]
as.numeric(res)
}
#Gr_x_k Gr_x_j g_i(x)
hefullg <- function(x, i, j, k)
{
x <- x[1:7]
if(i == 1)
{
he1 <- diag( c(0, 2, 6, 12, 20, 30, 42) * x ^ c(0, 0, 1:5) )
}
if(i == 2)
{
he1 <- matrix(0, 7, 7)
he1[1, -1] <- sapply(2:7, function(i) prod(x[-c(1, i)]))
he1[2, -2] <- sapply(c(1, 3:7), function(i) prod(x[-c(2, i)]))
he1[3, -3] <- sapply(c(1:2, 4:7), function(i) prod(x[-c(3, i)]))
he1[4, -4] <- sapply(c(1:3, 5:7), function(i) prod(x[-c(4, i)]))
he1[5, -5] <- sapply(c(1:4, 6:7), function(i) prod(x[-c(5, i)]))
he1[6, -6] <- sapply(c(1:5, 7:7), function(i) prod(x[-c(6, i)]))
he1[7, -7] <- sapply(1:6, function(i) prod(x[-c(7, i)]))
he2 <- matrix(0, 7, 7)
}
if(i == 3)
{
he1 <- diag(rep(2, 7))
he2 <- matrix(0, 7, 7)
}
if(i != 1)
return( c(he1[j, k], he2[j, k]) )
else
return( he1[j, k] )
}
# (3) compute Phi
#
z <- rexp(sum(dimx) + sum(dimlam))
n <- sum(dimx)
m <- sum(dimlam)
x <- z[1:n]
lam <- z[(n+1):(n+m)]
resphi <- funSSR(z, dimx, dimlam, grobj=grfullob, constr=g, grconstr=grfullg, compl=phiFB)
check <- c(grfullob(x, 1, 1) + lam[1] * grfullg(x, 1, 1),
grfullob(x, 1, 2) + lam[1] * grfullg(x, 1, 2),
grfullob(x, 2, 3) + lam[2:3] %*% grfullg(x, 2, 3),
grfullob(x, 2, 4) + lam[2:3] %*% grfullg(x, 2, 4),
grfullob(x, 3, 5) + lam[4:5] %*% grfullg(x, 3, 5),
grfullob(x, 3, 6) + lam[4:5] %*% grfullg(x, 3, 6),
grfullob(x, 3, 7) + lam[4:5] %*% grfullg(x, 3, 7),
phiFB( -g(x, 1), lam[1]),
phiFB( -g(x, 2)[1], lam[2]),
phiFB( -g(x, 2)[2], lam[3]),
phiFB( -g(x, 3)[1], lam[4]),
phiFB( -g(x, 3)[2], lam[5]))
#check
cat("\n\n________________________________________\n\n")
#part A
print(cbind(check, res=as.numeric(resphi))[1:n, ])
#part B
print(cbind(check, res=as.numeric(resphi))[(n+1):(n+m), ])
# (4) compute Jac Phi
#
resjacphi <- jacSSR(z, dimx, dimlam, heobj=hefullob, constr=g, grconstr=grfullg,
heconstr=hefullg, gcompla=GrAphiFB, gcomplb=GrBphiFB)
#check
cat("\n\n________________________________________\n\n")
cat("\n\npart A\n\n")
checkA <-
rbind(
c(hefullob(x, 1, 1, 1) + lam[1]*hefullg(x, 1, 1, 1),
hefullob(x, 1, 1, 2) + lam[1]*hefullg(x, 1, 1, 2),
hefullob(x, 1, 1, 3) + lam[1]*hefullg(x, 1, 1, 3),
hefullob(x, 1, 1, 4) + lam[1]*hefullg(x, 1, 1, 4),
hefullob(x, 1, 1, 5) + lam[1]*hefullg(x, 1, 1, 5),
hefullob(x, 1, 1, 6) + lam[1]*hefullg(x, 1, 1, 6),
hefullob(x, 1, 1, 7) + lam[1]*hefullg(x, 1, 1, 7)
),
c(hefullob(x, 1, 2, 1) + lam[1]*hefullg(x, 1, 2, 1),
hefullob(x, 1, 2, 2) + lam[1]*hefullg(x, 1, 2, 2),
hefullob(x, 1, 2, 3) + lam[1]*hefullg(x, 1, 2, 3),
hefullob(x, 1, 2, 4) + lam[1]*hefullg(x, 1, 2, 4),
hefullob(x, 1, 2, 5) + lam[1]*hefullg(x, 1, 2, 5),
hefullob(x, 1, 2, 6) + lam[1]*hefullg(x, 1, 2, 6),
hefullob(x, 1, 2, 7) + lam[1]*hefullg(x, 1, 2, 7)
),
c(hefullob(x, 2, 3, 1) + lam[2:3] %*% hefullg(x, 2, 3, 1),
hefullob(x, 2, 3, 2) + lam[2:3] %*% hefullg(x, 2, 3, 2),
hefullob(x, 2, 3, 3) + lam[2:3] %*% hefullg(x, 2, 3, 3),
hefullob(x, 2, 3, 4) + lam[2:3] %*% hefullg(x, 2, 3, 4),
hefullob(x, 2, 3, 5) + lam[2:3] %*% hefullg(x, 2, 3, 5),
hefullob(x, 2, 3, 6) + lam[2:3] %*% hefullg(x, 2, 3, 6),
hefullob(x, 2, 3, 7) + lam[2:3] %*% hefullg(x, 2, 3, 7)
),
c(hefullob(x, 2, 4, 1) + lam[2:3] %*% hefullg(x, 2, 4, 1),
hefullob(x, 2, 4, 2) + lam[2:3] %*% hefullg(x, 2, 4, 2),
hefullob(x, 2, 4, 3) + lam[2:3] %*% hefullg(x, 2, 4, 3),
hefullob(x, 2, 4, 4) + lam[2:3] %*% hefullg(x, 2, 4, 4),
hefullob(x, 2, 4, 5) + lam[2:3] %*% hefullg(x, 2, 4, 5),
hefullob(x, 2, 4, 6) + lam[2:3] %*% hefullg(x, 2, 4, 6),
hefullob(x, 2, 4, 7) + lam[2:3] %*% hefullg(x, 2, 4, 7)
),
c(hefullob(x, 3, 5, 1) + lam[4:5] %*% hefullg(x, 3, 5, 1),
hefullob(x, 3, 5, 2) + lam[4:5] %*% hefullg(x, 3, 5, 2),
hefullob(x, 3, 5, 3) + lam[4:5] %*% hefullg(x, 3, 5, 3),
hefullob(x, 3, 5, 4) + lam[4:5] %*% hefullg(x, 3, 5, 4),
hefullob(x, 3, 5, 5) + lam[4:5] %*% hefullg(x, 3, 5, 5),
hefullob(x, 3, 5, 6) + lam[4:5] %*% hefullg(x, 3, 5, 6),
hefullob(x, 3, 5, 7) + lam[4:5] %*% hefullg(x, 3, 5, 7)
),
c(hefullob(x, 3, 6, 1) + lam[4:5] %*% hefullg(x, 3, 6, 1),
hefullob(x, 3, 6, 2) + lam[4:5] %*% hefullg(x, 3, 6, 2),
hefullob(x, 3, 6, 3) + lam[4:5] %*% hefullg(x, 3, 6, 3),
hefullob(x, 3, 6, 4) + lam[4:5] %*% hefullg(x, 3, 6, 4),
hefullob(x, 3, 6, 5) + lam[4:5] %*% hefullg(x, 3, 6, 5),
hefullob(x, 3, 6, 6) + lam[4:5] %*% hefullg(x, 3, 6, 6),
hefullob(x, 3, 6, 7) + lam[4:5] %*% hefullg(x, 3, 6, 7)
),
c(hefullob(x, 3, 7, 1) + lam[4:5] %*% hefullg(x, 3, 7, 1),
hefullob(x, 3, 7, 2) + lam[4:5] %*% hefullg(x, 3, 7, 2),
hefullob(x, 3, 7, 3) + lam[4:5] %*% hefullg(x, 3, 7, 3),
hefullob(x, 3, 7, 4) + lam[4:5] %*% hefullg(x, 3, 7, 4),
hefullob(x, 3, 7, 5) + lam[4:5] %*% hefullg(x, 3, 7, 5),
hefullob(x, 3, 7, 6) + lam[4:5] %*% hefullg(x, 3, 7, 6),
hefullob(x, 3, 7, 7) + lam[4:5] %*% hefullg(x, 3, 7, 7)
)
)
print(resjacphi[1:n, 1:n] - checkA)
cat("\n\n________________________________________\n\n")
cat("\n\npart B\n\n")
checkB <-
rbind(
cbind(c(grfullg(x, 1, 1), grfullg(x, 1, 2)), c(0, 0), c(0, 0), c(0, 0), c(0, 0)),
cbind(c(0, 0), rbind(grfullg(x, 2, 3), grfullg(x, 2, 4)), c(0, 0), c(0, 0)),
cbind(c(0, 0, 0), c(0, 0, 0), c(0, 0, 0),
rbind(grfullg(x, 3, 5), grfullg(x, 3, 6), grfullg(x, 3, 7)))
)
print(resjacphi[1:n, (n+1):(n+m)] - checkB)
cat("\n\n________________________________________\n\n")
cat("\n\npart C\n\n")
gx <- c(g(x,1), g(x,2), g(x,3))
checkC <-
- t(
cbind(
rbind(
grfullg(x, 1, 1) * GrAphiFB(-gx, lam)[1],
grfullg(x, 1, 2) * GrAphiFB(-gx, lam)[1],
grfullg(x, 1, 3) * GrAphiFB(-gx, lam)[1],
grfullg(x, 1, 4) * GrAphiFB(-gx, lam)[1],
grfullg(x, 1, 5) * GrAphiFB(-gx, lam)[1],
grfullg(x, 1, 6) * GrAphiFB(-gx, lam)[1],
grfullg(x, 1, 7) * GrAphiFB(-gx, lam)[1]
),
rbind(
grfullg(x, 2, 1) * GrAphiFB(-gx, lam)[2:3],
grfullg(x, 2, 2) * GrAphiFB(-gx, lam)[2:3],
grfullg(x, 2, 3) * GrAphiFB(-gx, lam)[2:3],
grfullg(x, 2, 4) * GrAphiFB(-gx, lam)[2:3],
grfullg(x, 2, 5) * GrAphiFB(-gx, lam)[2:3],
grfullg(x, 2, 6) * GrAphiFB(-gx, lam)[2:3],
grfullg(x, 2, 7) * GrAphiFB(-gx, lam)[2:3]
),
rbind(
grfullg(x, 3, 1) * GrAphiFB(-gx, lam)[4:5],
grfullg(x, 3, 2) * GrAphiFB(-gx, lam)[4:5],
grfullg(x, 3, 3) * GrAphiFB(-gx, lam)[4:5],
grfullg(x, 3, 4) * GrAphiFB(-gx, lam)[4:5],
grfullg(x, 3, 5) * GrAphiFB(-gx, lam)[4:5],
grfullg(x, 3, 6) * GrAphiFB(-gx, lam)[4:5],
grfullg(x, 3, 7) * GrAphiFB(-gx, lam)[4:5]
)
)
)
print(resjacphi[(n+1):(n+m), 1:n] - checkC)
cat("\n\n________________________________________\n\n")
cat("\n\npart D\n\n")
checkD <- diag(GrBphiFB(-gx, lam))
print(resjacphi[(n+1):(n+m), (n+1):(n+m)] - checkD)
GNE documentation built on Oct. 17, 2024, 5:08 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| SSR | R Documentation |
SemiSmooth Reformulation
Description
functions of the SemiSmooth Reformulation of the GNEP
Usage
funSSR(z, dimx, dimlam, grobj, arggrobj, constr, argconstr, grconstr, arggrconstr,
compl, argcompl, dimmu, joint, argjoint, grjoint, arggrjoint, echo=FALSE)
jacSSR(z, dimx, dimlam, heobj, argheobj, constr, argconstr, grconstr, arggrconstr,
heconstr, argheconstr, gcompla, gcomplb, argcompl, dimmu, joint, argjoint,
grjoint, arggrjoint, hejoint, arghejoint, echo=FALSE)
Arguments
z |
a numeric vector |
dimx |
a vector of dimension for |
dimlam |
a vector of dimension for |
grobj |
gradient of the objective function, see details. |
arggrobj |
a list of additional arguments of the objective gradient. |
constr |
constraint function, see details. |
argconstr |
a list of additional arguments of the constraint function. |
grconstr |
gradient of the constraint function, see details. |
arggrconstr |
a list of additional arguments of the constraint gradient. |
compl |
the complementarity function with (at least) two arguments: |
argcompl |
list of possible additional arguments for |
dimmu |
a vector of dimension for |
joint |
joint function, see details. |
argjoint |
a list of additional arguments of the joint function. |
grjoint |
gradient of the joint function, see details. |
arggrjoint |
a list of additional arguments of the joint gradient. |
heobj |
Hessian of the objective function, see details. |
argheobj |
a list of additional arguments of the objective Hessian. |
heconstr |
Hessian of the constraint function, see details. |
argheconstr |
a list of additional arguments of the constraint Hessian. |
gcompla |
derivative of the complementarity function w.r.t. the first argument. |
gcomplb |
derivative of the complementarity function w.r.t. the second argument. |
hejoint |
Hessian of the joint function, see details. |
arghejoint |
a list of additional arguments of the joint Hessian. |
echo |
a logical to show some traces. |
Details
Compute the SemiSmooth Reformulation of the GNEP: the Generalized Nash equilibrium problem is defined
by objective functions Obj with player variables x defined in dimx and
may have player-dependent constraint functions g of dimension dimlam
and/or a common shared joint function h of dimension dimmu,
where the Lagrange multiplier are lambda and mu, respectively,
see F. Facchinei et al.(2009) where there is no joint function.
- Arguments of the Phi function
-
The arguments which are functions must respect the following features
grobj-
The gradient
Grad Objof an objective functionObj(to be minimized) must have 3 arguments forGrad Obj(z, playnum, ideriv): vectorz, player number, derivative index , and optionnally additional arguments inarggrobj. constr-
The constraint function
gmust have 2 arguments: vectorz, player number, such thatg(z, playnum) <= 0. Optionnally,gmay have additional arguments inargconstr. grconstr-
The gradient of the constraint function
gmust have 3 arguments: vectorz, player number, derivative index, and optionnally additional arguments inarggrconstr. complIt must have two arguments and optionnally additional arguments in
argcompl. A typical example is the minimum function.joint-
The constraint function
hmust have 1 argument: vectorz, such thath(z) <= 0. Optionnally,hmay have additional arguments inargjoint. grjoint-
The gradient of the constraint function
hmust have 2 arguments: vectorz, derivative index, and optionnally additional arguments inarggrjoint.
- Arguments of the Jacobian of Phi
-
The arguments which are functions must respect the following features
heobjIt must have 4 arguments: vector
z, player number, two derivative indexes and optionnally additional arguments inargheobj.heconstrIt must have 4 arguments: vector
z, player number, two derivative indexes and optionnally additional arguments inargheconstr.gcompla,gcomplbIt must have two arguments and optionnally additional arguments in
argcompl.hejointIt must have 3 arguments: vector
z, two derivative indexes and optionnally additional arguments inarghejoint.
See the example below.
Value
A vector for funSSR or a matrix for jacSSR.
Author(s)
Christophe Dutang
References
F. Facchinei, A. Fischer and V. Piccialli (2009), Generalized Nash equilibrium problems and Newton methods, Math. Program.
See Also
See also GNE.nseq.
Examples
# (1) associated objective functions
#
dimx <- c(2, 2, 3)
#Gr_x_j O_i(x)
grfullob <- function(x, i, j)
{
x <- x[1:7]
if(i == 1)
{
grad <- 3*(x - 1:7)^2
}
if(i == 2)
{
grad <- 1:7*(x - 1:7)^(0:6)
}
if(i == 3)
{
s <- x[5]^2 + x[6]^2 + x[7]^2 - 5
grad <- c(1, 0, 1, 0, 4*x[5]*s, 4*x[6]*s, 4*x[7]*s)
}
grad[j]
}
#Gr_x_k Gr_x_j O_i(x)
hefullob <- function(x, i, j, k)
{
x <- x[1:7]
if(i == 1)
{
he <- diag( 6*(x - 1:7) )
}
if(i == 2)
{
he <- diag( c(0, 2, 6, 12, 20, 30, 42)*(x - 1:7)^c(0, 0:5) )
}
if(i == 3)
{
s <- x[5]^2 + x[6]^2 + x[7]^2
he <- rbind(rep(0, 7), rep(0, 7), rep(0, 7), rep(0, 7),
c(0, 0, 0, 0, 4*s+8*x[5]^2, 8*x[5]*x[6], 8*x[5]*x[7]),
c(0, 0, 0, 0, 8*x[5]*x[6], 4*s+8*x[6]^2, 8*x[6]*x[7]),
c(0, 0, 0, 0, 8*x[5]*x[7], 8*x[6]*x[7], 4*s+8*x[7]^2) )
}
he[j,k]
}
# (2) constraint linked functions
#
dimlam <- c(1, 2, 2)
#constraint function g_i(x)
g <- function(x, i)
{
x <- x[1:7]
if(i == 1)
res <- sum( x^(1:7) ) -7
if(i == 2)
res <- c(sum(x) + prod(x) - 14, 20 - sum(x))
if(i == 3)
res <- c(sum(x^2) + 1, 100 - sum(x))
res
}
#Gr_x_j g_i(x)
grfullg <- function(x, i, j)
{
x <- x[1:7]
if(i == 1)
{
grad <- (1:7) * x ^ (0:6)
}
if(i == 2)
{
grad <- 1 + sapply(1:7, function(i) prod(x[-i]))
grad <- cbind(grad, -1)
}
if(i == 3)
{
grad <- cbind(2*x, -1)
}
if(i == 1)
res <- grad[j]
if(i != 1)
res <- grad[j,]
as.numeric(res)
}
#Gr_x_k Gr_x_j g_i(x)
hefullg <- function(x, i, j, k)
{
x <- x[1:7]
if(i == 1)
{
he1 <- diag( c(0, 2, 6, 12, 20, 30, 42) * x ^ c(0, 0, 1:5) )
}
if(i == 2)
{
he1 <- matrix(0, 7, 7)
he1[1, -1] <- sapply(2:7, function(i) prod(x[-c(1, i)]))
he1[2, -2] <- sapply(c(1, 3:7), function(i) prod(x[-c(2, i)]))
he1[3, -3] <- sapply(c(1:2, 4:7), function(i) prod(x[-c(3, i)]))
he1[4, -4] <- sapply(c(1:3, 5:7), function(i) prod(x[-c(4, i)]))
he1[5, -5] <- sapply(c(1:4, 6:7), function(i) prod(x[-c(5, i)]))
he1[6, -6] <- sapply(c(1:5, 7:7), function(i) prod(x[-c(6, i)]))
he1[7, -7] <- sapply(1:6, function(i) prod(x[-c(7, i)]))
he2 <- matrix(0, 7, 7)
}
if(i == 3)
{
he1 <- diag(rep(2, 7))
he2 <- matrix(0, 7, 7)
}
if(i != 1)
return( c(he1[j, k], he2[j, k]) )
else
return( he1[j, k] )
}
# (3) compute Phi
#
z <- rexp(sum(dimx) + sum(dimlam))
n <- sum(dimx)
m <- sum(dimlam)
x <- z[1:n]
lam <- z[(n+1):(n+m)]
resphi <- funSSR(z, dimx, dimlam, grobj=grfullob, constr=g, grconstr=grfullg, compl=phiFB)
check <- c(grfullob(x, 1, 1) + lam[1] * grfullg(x, 1, 1),
grfullob(x, 1, 2) + lam[1] * grfullg(x, 1, 2),
grfullob(x, 2, 3) + lam[2:3] %*% grfullg(x, 2, 3),
grfullob(x, 2, 4) + lam[2:3] %*% grfullg(x, 2, 4),
grfullob(x, 3, 5) + lam[4:5] %*% grfullg(x, 3, 5),
grfullob(x, 3, 6) + lam[4:5] %*% grfullg(x, 3, 6),
grfullob(x, 3, 7) + lam[4:5] %*% grfullg(x, 3, 7),
phiFB( -g(x, 1), lam[1]),
phiFB( -g(x, 2)[1], lam[2]),
phiFB( -g(x, 2)[2], lam[3]),
phiFB( -g(x, 3)[1], lam[4]),
phiFB( -g(x, 3)[2], lam[5]))
#check
cat("\n\n________________________________________\n\n")
#part A
print(cbind(check, res=as.numeric(resphi))[1:n, ])
#part B
print(cbind(check, res=as.numeric(resphi))[(n+1):(n+m), ])
# (4) compute Jac Phi
#
resjacphi <- jacSSR(z, dimx, dimlam, heobj=hefullob, constr=g, grconstr=grfullg,
heconstr=hefullg, gcompla=GrAphiFB, gcomplb=GrBphiFB)
#check
cat("\n\n________________________________________\n\n")
cat("\n\npart A\n\n")
checkA <-
rbind(
c(hefullob(x, 1, 1, 1) + lam[1]*hefullg(x, 1, 1, 1),
hefullob(x, 1, 1, 2) + lam[1]*hefullg(x, 1, 1, 2),
hefullob(x, 1, 1, 3) + lam[1]*hefullg(x, 1, 1, 3),
hefullob(x, 1, 1, 4) + lam[1]*hefullg(x, 1, 1, 4),
hefullob(x, 1, 1, 5) + lam[1]*hefullg(x, 1, 1, 5),
hefullob(x, 1, 1, 6) + lam[1]*hefullg(x, 1, 1, 6),
hefullob(x, 1, 1, 7) + lam[1]*hefullg(x, 1, 1, 7)
),
c(hefullob(x, 1, 2, 1) + lam[1]*hefullg(x, 1, 2, 1),
hefullob(x, 1, 2, 2) + lam[1]*hefullg(x, 1, 2, 2),
hefullob(x, 1, 2, 3) + lam[1]*hefullg(x, 1, 2, 3),
hefullob(x, 1, 2, 4) + lam[1]*hefullg(x, 1, 2, 4),
hefullob(x, 1, 2, 5) + lam[1]*hefullg(x, 1, 2, 5),
hefullob(x, 1, 2, 6) + lam[1]*hefullg(x, 1, 2, 6),
hefullob(x, 1, 2, 7) + lam[1]*hefullg(x, 1, 2, 7)
),
c(hefullob(x, 2, 3, 1) + lam[2:3] %*% hefullg(x, 2, 3, 1),
hefullob(x, 2, 3, 2) + lam[2:3] %*% hefullg(x, 2, 3, 2),
hefullob(x, 2, 3, 3) + lam[2:3] %*% hefullg(x, 2, 3, 3),
hefullob(x, 2, 3, 4) + lam[2:3] %*% hefullg(x, 2, 3, 4),
hefullob(x, 2, 3, 5) + lam[2:3] %*% hefullg(x, 2, 3, 5),
hefullob(x, 2, 3, 6) + lam[2:3] %*% hefullg(x, 2, 3, 6),
hefullob(x, 2, 3, 7) + lam[2:3] %*% hefullg(x, 2, 3, 7)
),
c(hefullob(x, 2, 4, 1) + lam[2:3] %*% hefullg(x, 2, 4, 1),
hefullob(x, 2, 4, 2) + lam[2:3] %*% hefullg(x, 2, 4, 2),
hefullob(x, 2, 4, 3) + lam[2:3] %*% hefullg(x, 2, 4, 3),
hefullob(x, 2, 4, 4) + lam[2:3] %*% hefullg(x, 2, 4, 4),
hefullob(x, 2, 4, 5) + lam[2:3] %*% hefullg(x, 2, 4, 5),
hefullob(x, 2, 4, 6) + lam[2:3] %*% hefullg(x, 2, 4, 6),
hefullob(x, 2, 4, 7) + lam[2:3] %*% hefullg(x, 2, 4, 7)
),
c(hefullob(x, 3, 5, 1) + lam[4:5] %*% hefullg(x, 3, 5, 1),
hefullob(x, 3, 5, 2) + lam[4:5] %*% hefullg(x, 3, 5, 2),
hefullob(x, 3, 5, 3) + lam[4:5] %*% hefullg(x, 3, 5, 3),
hefullob(x, 3, 5, 4) + lam[4:5] %*% hefullg(x, 3, 5, 4),
hefullob(x, 3, 5, 5) + lam[4:5] %*% hefullg(x, 3, 5, 5),
hefullob(x, 3, 5, 6) + lam[4:5] %*% hefullg(x, 3, 5, 6),
hefullob(x, 3, 5, 7) + lam[4:5] %*% hefullg(x, 3, 5, 7)
),
c(hefullob(x, 3, 6, 1) + lam[4:5] %*% hefullg(x, 3, 6, 1),
hefullob(x, 3, 6, 2) + lam[4:5] %*% hefullg(x, 3, 6, 2),
hefullob(x, 3, 6, 3) + lam[4:5] %*% hefullg(x, 3, 6, 3),
hefullob(x, 3, 6, 4) + lam[4:5] %*% hefullg(x, 3, 6, 4),
hefullob(x, 3, 6, 5) + lam[4:5] %*% hefullg(x, 3, 6, 5),
hefullob(x, 3, 6, 6) + lam[4:5] %*% hefullg(x, 3, 6, 6),
hefullob(x, 3, 6, 7) + lam[4:5] %*% hefullg(x, 3, 6, 7)
),
c(hefullob(x, 3, 7, 1) + lam[4:5] %*% hefullg(x, 3, 7, 1),
hefullob(x, 3, 7, 2) + lam[4:5] %*% hefullg(x, 3, 7, 2),
hefullob(x, 3, 7, 3) + lam[4:5] %*% hefullg(x, 3, 7, 3),
hefullob(x, 3, 7, 4) + lam[4:5] %*% hefullg(x, 3, 7, 4),
hefullob(x, 3, 7, 5) + lam[4:5] %*% hefullg(x, 3, 7, 5),
hefullob(x, 3, 7, 6) + lam[4:5] %*% hefullg(x, 3, 7, 6),
hefullob(x, 3, 7, 7) + lam[4:5] %*% hefullg(x, 3, 7, 7)
)
)
print(resjacphi[1:n, 1:n] - checkA)
cat("\n\n________________________________________\n\n")
cat("\n\npart B\n\n")
checkB <-
rbind(
cbind(c(grfullg(x, 1, 1), grfullg(x, 1, 2)), c(0, 0), c(0, 0), c(0, 0), c(0, 0)),
cbind(c(0, 0), rbind(grfullg(x, 2, 3), grfullg(x, 2, 4)), c(0, 0), c(0, 0)),
cbind(c(0, 0, 0), c(0, 0, 0), c(0, 0, 0),
rbind(grfullg(x, 3, 5), grfullg(x, 3, 6), grfullg(x, 3, 7)))
)
print(resjacphi[1:n, (n+1):(n+m)] - checkB)
cat("\n\n________________________________________\n\n")
cat("\n\npart C\n\n")
gx <- c(g(x,1), g(x,2), g(x,3))
checkC <-
- t(
cbind(
rbind(
grfullg(x, 1, 1) * GrAphiFB(-gx, lam)[1],
grfullg(x, 1, 2) * GrAphiFB(-gx, lam)[1],
grfullg(x, 1, 3) * GrAphiFB(-gx, lam)[1],
grfullg(x, 1, 4) * GrAphiFB(-gx, lam)[1],
grfullg(x, 1, 5) * GrAphiFB(-gx, lam)[1],
grfullg(x, 1, 6) * GrAphiFB(-gx, lam)[1],
grfullg(x, 1, 7) * GrAphiFB(-gx, lam)[1]
),
rbind(
grfullg(x, 2, 1) * GrAphiFB(-gx, lam)[2:3],
grfullg(x, 2, 2) * GrAphiFB(-gx, lam)[2:3],
grfullg(x, 2, 3) * GrAphiFB(-gx, lam)[2:3],
grfullg(x, 2, 4) * GrAphiFB(-gx, lam)[2:3],
grfullg(x, 2, 5) * GrAphiFB(-gx, lam)[2:3],
grfullg(x, 2, 6) * GrAphiFB(-gx, lam)[2:3],
grfullg(x, 2, 7) * GrAphiFB(-gx, lam)[2:3]
),
rbind(
grfullg(x, 3, 1) * GrAphiFB(-gx, lam)[4:5],
grfullg(x, 3, 2) * GrAphiFB(-gx, lam)[4:5],
grfullg(x, 3, 3) * GrAphiFB(-gx, lam)[4:5],
grfullg(x, 3, 4) * GrAphiFB(-gx, lam)[4:5],
grfullg(x, 3, 5) * GrAphiFB(-gx, lam)[4:5],
grfullg(x, 3, 6) * GrAphiFB(-gx, lam)[4:5],
grfullg(x, 3, 7) * GrAphiFB(-gx, lam)[4:5]
)
)
)
print(resjacphi[(n+1):(n+m), 1:n] - checkC)
cat("\n\n________________________________________\n\n")
cat("\n\npart D\n\n")
checkD <- diag(GrBphiFB(-gx, lam))
print(resjacphi[(n+1):(n+m), (n+1):(n+m)] - checkD)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.