Last data update: 2014.03.03

 R: SemiSmooth Reformulation
SSRR 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.

compl

It 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

heobj

It must have 4 arguments: vector z, player number, two derivative indexes and optionnally additional arguments in argheobj.

heconstr

It must have 4 arguments: vector z, player number, two derivative indexes and optionnally additional arguments in argheconstr.

gcompla,gcomplb

It must have two arguments and optionnally additional arguments in argcompl.

hejoint

It 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)

Results


R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(GNE)
Loading required package: alabama
Loading required package: numDeriv
Loading required package: nleqslv
Loading required package: BB
Loading required package: SQUAREM
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/GNE/util-SSR.Rd_%03d_medium.png", width=480, height=480)
> ### Name: SSR
> ### Title: SemiSmooth Reformulation
> ### Aliases: SSR funSSR jacSSR
> ### Keywords: math optimize
> 
> ### ** 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, ])
            check          res
[1,]    4.7335070    4.7335070
[2,]    8.9324574    8.9324574
[3,]   12.9928537   12.9928537
[4,] -257.1265856 -257.1265856
[5,]   -0.1500646   -0.1500646
[6,]   18.9785509   18.9785509
[7,]   -2.5117080   -2.5117080
> #part B
> print(cbind(check, res=as.numeric(resphi))[(n+1):(n+m), ])
           check         res
[1,] 649.4925817 649.4925817
[2,]  -0.5154893  -0.5154893
[3,]  27.6247860  27.6247860
[4,]  18.0286559  18.0286559
[5,] 187.6865607 187.6865607
> 	
> # (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")	


part A

> 
> 
> 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)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    0    0    0    0    0    0    0
[2,]    0    0    0    0    0    0    0
[3,]    0    0    0    0    0    0    0
[4,]    0    0    0    0    0    0    0
[5,]    0    0    0    0    0    0    0
[6,]    0    0    0    0    0    0    0
[7,]    0    0    0    0    0    0    0
> 
> 
> cat("\n\n________________________________________\n\n")


________________________________________

> 
> 
> cat("\n\npart B\n\n")	


part B

> 
> 
> 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)	
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    0    0    0    0
[2,]    0    0    0    0    0
[3,]    0    0    0    0    0
[4,]    0    0    0    0    0
[5,]    0    0    0    0    0
[6,]    0    0    0    0    0
[7,]    0    0    0    0    0
> 
> 
> cat("\n\n________________________________________\n\n")


________________________________________

> cat("\n\npart C\n\n")	


part C

> 
> 
> 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)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    0    0    0    0    0    0    0
[2,]    0    0    0    0    0    0    0
[3,]    0    0    0    0    0    0    0
[4,]    0    0    0    0    0    0    0
[5,]    0    0    0    0    0    0    0
> 
> 
> cat("\n\n________________________________________\n\n")


________________________________________

> 
> cat("\n\npart D\n\n")	


part D

> 
> 
> checkD <- diag(GrBphiFB(-gx, lam)) 
> 
> print(resjacphi[(n+1):(n+m), (n+1):(n+m)] - checkD)
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    0    0    0    0
[2,]    0    0    0    0    0
[3,]    0    0    0    0    0
[4,]    0    0    0    0    0
[5,]    0    0    0    0    0
> 
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>


Created & Maintained by Osamu Ogasawara (osamu.ogasawara@gmail.com) and IMSBIO