Nothing
#########################################################
### Authors: Giulia Marcon and Simone Padoan ###
### Emails: giulia.marcon@phd.unibocconi.it, ###
### simone.padoan@unibocconi.it ###
### Institution: Department of Decision Sciences, ###
### University Bocconi of Milan ###
### File name: ConvexityConstraints.r ###
### Description: ###
### This file contains a set of procedures ###
### for computing the matrix of convexity constraints ###
### used in the quadprog procedure of optimization ###
### Last change: 15/08/2016 ###
#########################################################
## Permutation function taken from gtools package (hidden function)
# Main routine that provides a matrix of CONSTRAINTS that garanties
# the CONVEXITY.
convexity <-function(d, k){
if(d==2){
R1 <- matrix(0,k-1,k+1)
for (i in 1:(k-1))
R1[i,i:(i+2)] <- c(1,-2,1)
}
else{
###########################
#Subroutines for convexity:
# 1)
index <- function(k, d){
beta.index <- expand.grid(rep(list(0:k),d))
beta.index.sum <- rowSums(beta.index)
restr <- which(beta.index.sum<=k)
v <- as.matrix(beta.index[restr,ncol(beta.index):1])
return(v)
}
# 2)
lookUpTable <- function(kTilde, dTilde){
nrows= (kTilde+1)^(dTilde)
lookup <- matrix(nrow = nrows,ncol=2)
feasibles =0;
for(i in 0:nrows-1){
if(isFeasible(i,kTilde+1)){
lookup[i+1,1] = i
lookup[i+1,2] = feasibles
feasibles = feasibles+1;
}
}
return(lookup)
}
# 3)
isFeasible <- function (number, base){
dividend = number
sum =0
while(dividend>=base){
sum = sum + dividend%%base
dividend = floor(dividend/base)
}
sum =sum+dividend
return(sum<=(base-1))
}
# 4)
frombaseKToBase10 <-function(row, kTilde){
dTilde <- length(row)
baseNum <- kTilde+1
index <- 0;
for(i in dTilde:1){
index = index + row[i]*baseNum^(dTilde-i)
}
return(index)
}
# 5)
numberToModifyInCompleteTable <- function(numberBase10, power_direction, kTilde, order){
return(numberBase10 + ((kTilde+1)^(power_direction-1))*order)
}
# 6)
diff1inDirections <-function (lookUp,numberBase10, k, vLength, power_direction1, power_direction2) {
resultVector = vector(length = vLength,mode = "numeric")
numberAfterFirstDirection = numberToModifyInCompleteTable(numberBase10, power_direction1, k,1)
indexToBeModified11Complete = numberToModifyInCompleteTable(numberAfterFirstDirection, power_direction2, k,1) + 1
index11 = lookUp[indexToBeModified11Complete,2]
indexToBeModified01Complete = numberToModifyInCompleteTable(numberBase10, power_direction2, k,1) + 1
index01 = lookUp[indexToBeModified01Complete,2]
indexToBeModified10Complete = numberToModifyInCompleteTable(numberBase10, power_direction1, k,1) + 1
index10 = lookUp[indexToBeModified10Complete,2]
index0 = lookUp[numberBase10 +1,2]
resultVector[index11] = 1;
resultVector[index10] = resultVector[index10]-1;
resultVector[index01] = resultVector[index01] -1;
resultVector[index0] = 1;
return(resultVector)
}
##########################################
dTilde = d-1
kTilde = k-2
v = index(k, dTilde)
vSub = index(kTilde, dTilde)
vLengthSub = dim(vSub)[1]
vLength = dim(v)[1]
lookUp = lookUpTable(k, dTilde)
directions = 1: (dTilde)
perm <- permutations(2,dTilde-1,c(1,-1),repeats.allowed =TRUE);
dimPerm <- dim(perm)[1]
R1 = matrix(data= 0, nrow = dTilde*dimPerm*vLengthSub, ncol=vLength)
for(vIndex in 1:vLengthSub){
numberBase10 = frombaseKToBase10(vSub[vIndex,],k)
for(i in directions){
diff1is = matrix(nrow = dTilde, ncol=vLength)
diff2s = matrix(nrow = dTilde, ncol=vLength)
diff2s = diff1inDirections(lookUp,numberBase10, k, vLength, i, i)
spuriousDirections = directions[-i]
for(j in spuriousDirections){
diff1is[j,] =diff1inDirections(lookUp,numberBase10, k, vLength, i, j)
}
for(t in 1:dimPerm){
begin =(vIndex-1)*dimPerm*dTilde
rowR1 = diff2s + perm[t,] %*%diff1is[-i,]
R1[begin +(dTilde-i)*dimPerm+t,] = rowR1
}
}
}
}
return(R1)
}
# Main routine that provides a matrix of CONSTRAINTS that garanties
# the CONVEXITY. OLD VERSION
convexity_old <- function(v, d){
combn_eqk <- v
q <- nrow(v)
k <- max(v)
if(d==2){
Aconvx <- matrix(0,k-1,k+1)
for (i in 1:(k-1))
Aconvx[i,i:(i+2)] <- c(1,-2,1)
}
if(d==3){
B <- matrix(c(0, 1, rep(0,k-1), -1, -1, rep(0,k-1), 1,
2, -1, rep(0,k-1), -3, 1, rep(0,k-1), 1,
0, -1, 1, rep(0,k-2), 1, -1, rep(0,k),
2, -3, 1, rep(0,k-2), -1, 1, rep(0,k)),
4,ncol=3+2*k,byrow=T)
ri <- 0
co <- 0
for (i in 1:(k-1)){
ri <- ri + k-i
co <- c(co, (i+(i-1)*k):(i*k-1))
}
righe <- seq(1,ri*4,by=4)
colonne <- co[-1]
Aconvx <- matrix(0,length(righe)*nrow(B),(k+1)^2)
for (i in 1:length(righe)){
Aconvx[righe[i]:(righe[i]+nrow(B)-1),colonne[i]:(colonne[i]+ncol(B)-1)] <- B
}
restr <- labels(v)[[1]]
restr <- as.numeric(restr)
Aconvx <- Aconvx[,restr]
}
if(d>3){
sumvb <- apply(v,1,sum)
leqk_2 <- which(sumvb<=k-2)
combn_leqk_2 <- as.matrix(v[leqk_2,])
##########################
# Subroutines of Convexity:
# 1) Difference operators for convexity constraints
DELTAl0 <- function(d,combn_leqk_2,diffl,diff0,l){
diffll <- diffl0 <- diffl
diff0l <- diff00 <- diff0
for(m in 1:(d-1)){
for (i in 1:nrow(combn_leqk_2)){
diffll[i,l+1,m] <- diffl[i,l+1,m] + 1
diffl0[i,1,m] <- diffl[i,1,m] + 1
diff0l[i,l+1,m] <- diff0[i,l+1,m] + 1
diff00[i,1,m] <- diff0[i,1,m] + 1
}
}
return(list(diffll=diffll,diffl0=diffl0,
diff0l=diff0l,diff00=diff00))
}
# 2)
DELTAl0d <- function(d,combn_leqk_2,diffl,diff0,l){
diffll <- diffl0 <- diffl
diff0l <- diff00 <- diff0
for(m in 1:(d-1)){
for (i in 1:nrow(combn_leqk_2)){
diffll[i,l,m] <- diffl[i,l,m] + 1
diffl0[i,1,m] <- diffl[i,1,m]
diff0l[i,l,m] <- diff0[i,l,m] + 1
diff00[i,1,m] <- diff0[i,1,m]
}
}
return(list(diffll=diffll,diffl0=diffl0,
diff0l=diff0l,diff00=diff00))
}
# 3)
DELTAl0DELTAl0 <- function(d,combn_leqk_2,diffl,diff0,l){
deltal0 <- NULL
for (l in 1:(d-1)){
deltal0[[l]] <- DELTAl0(d=d,combn_leqk_2,diffl,diff0,l=l)
}
return(deltal0)
}
# 4)
DELTAl0DELTAl0d <- function(d,combn_leqk_2,diffl,diff0,l){
deltal0 <- NULL
for (l in 1:(d-1)){
deltal0[[l]] <- DELTAl0d(d=d,combn_leqk_2,diffl,diff0,l=l)
}
return(deltal0)
}
# DELTA_i0
diffl <- diff0 <- array(rep(combn_leqk_2,length(leqk_2)),dim=c(length(leqk_2),d-1,d-1))
for (l in 1:(d-1)){
for (i in 1:length(leqk_2)){
diffl[i,l,l] <- combn_leqk_2[i,l] + 1
diff0[i,1,l] <- combn_leqk_2[i,1]
}
}
# DELTAl0 DELTAl0
DIFF <- DELTAl0DELTAl0d(d=d,combn_leqk_2,diffl,diff0,l=l)
# + and - combinations
perm <- permutations(2,(d-2),c('+','-'),repeats.allowed=TRUE)
piuomeno <- NULL
piuomeno[[1]] <- cbind(rep("+",nrow(perm)),perm)
piuomeno[[d-1]] <- cbind(perm,rep("+",nrow(perm)))
for(l in 1:(d-3)){
piuomeno[[l+1]] <- cbind(perm[,1:l],rep("+",nrow(perm)),perm[,(l+1):(d-2)])
}
rowA <- nrow(perm)*nrow(combn_leqk_2)*(d-1)
Aconvx <- matrix(0,rowA,q)
ind <-0
f <- 0
o <- 0
g <- 0
for (l in 1:(d-1)){
g <- (l-1) * nrow(perm)
for(m in 1:(d-1)){
for(i in 1:nrow(combn_leqk_2)){
f <- (i-1+o) * nrow(perm)
for(h in 1:(nrow(perm))){
for(j in 1:nrow(combn_eqk)){
ind <- h+f+g
if(l==m){
if(all(DIFF[[l]]$diffll[i,,m]==combn_eqk[j,]))
Aconvx[ind,j] <- Aconvx[ind,j] + 1
if(all(DIFF[[l]]$diffl0[i,,m]==combn_eqk[j,]))
Aconvx[ind,j] <- Aconvx[ind,j] - 1
if(all(DIFF[[l]]$diff0l[i,,m]==combn_eqk[j,]))
Aconvx[ind,j] <- Aconvx[ind,j] - 1
if(all(DIFF[[l]]$diff00[i,,m]==combn_eqk[j,]))
Aconvx[ind,j] <- Aconvx[ind,j] + 1
}
else{
pm <- switch(piuomeno[[m]][h,l],"+"= 1, "-"= -1)
if(all(DIFF[[l]]$diffll[i,,m]==combn_eqk[j,]))
Aconvx[ind,j] <- Aconvx[ind,j] + 1 * pm
if(all(DIFF[[l]]$diffl0[i,,m]==combn_eqk[j,]))
Aconvx[ind,j] <- Aconvx[ind,j] - 1 * pm
if(all(DIFF[[l]]$diff0l[i,,m]==combn_eqk[j,]))
Aconvx[ind,j] <- Aconvx[ind,j] - 1 * pm
if(all(DIFF[[l]]$diff00[i,,m]==combn_eqk[j,]))
Aconvx[ind,j] <- Aconvx[ind,j] + 1 * pm
}
}
}
}
}
o <- o + 2
}
}
return(Aconvx)
}
###################################################################
########################### START CONSTRAINTS #####################
########################### ONLY BIVARIATE CASE ###################
constraints <- function(k){
kp <- k+1
km <- k-1
# Create the A, b0 matrix for constraints
# R1) Convexity
R1 <- matrix(0,km,kp)
for (i in 1:km)
R1[i,i:(i+2)] <- c(1,-2,1)
r1 <- rep(0,nrow(R1))
# R2) upper bound
# A(ei) = 1
R2 <- matrix(0,4,kp)
R2[1,1] <- R2[2,kp] <- 1
R2[3,1]<- R2[4,kp] <- -1
r2 <- c(1,1,-1,-1)
# R3) lower bound
# A(w) >= max(w)
R3 <- matrix(0,2,kp)
R3[1,2] <- R3[2,k] <- 1
r3 <- rep(1-1/k,2)
return(list(R=rbind(R1, R2, R3), r = c(r1, r2, r3) ))
}
########################### ONLY BIVARIATE CASE ###################
########################## END CONSTRAINTS ########################
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.