[S] Q on an optimization program (Simplex Method)

Mikio Kaihara (mkaihara@ichinoseki.ac.jp)
Mon, 29 Jun 1998 08:22:33 +0900


Dear Sir,
I(a beginner from May, Win 95, S-PLUS ver.4.0) am now fabricating a function,
"fok" which will minimize the "uu" by "Simplex Method".
"uu" was calculated by the function, "fun". The number of variables is 4, which
are the factors of the matrix (2X2), "tt(=ab[,,1:5])".
The function, "fun1" determine the 5 default values of "tt (=ab[,,1:5])" and
calculated 5 "uu". (We have 4 variables in "tt"
and need 5(4+1) points in case of "Simplex Method")
The function, "fun2" order the "tt(=ab[,,1:5])" and the "uu" in the descending
order.
The function, "fun3" change the first "tt(=ab[,,1])", which gives the largest
"uu(=u[1,1])" and calculate the new "uu(=new u[1,1]).
The first "tt(=ab[,,1])" will be calculated from the second to fifth
"tt(=ab[,,2---5])".
I changed the iteration number "S" of "fok", but the results were
circulating(???). And I was not able to find the reason.
Would you please help me ?
----------------------------------------------------------------------
> fok(20,9,2,1)
> uu
[,1]
[1,] 5.566436
[2,] 1718.825011
[3,] 1302.940425
[4,] 751.572541
[5,] 230.956874
> ab

, , 1
[,1] [,2]
[1,] 0.7 -0.3
[2,] -0.3 0.7

, , 2
[,1] [,2]
[1,] 1.9 0.9
[2,] 0.9 1.9

, , 3
[,1] [,2]
[1,] 1.6 0.6
[2,] 0.6 1.6

, , 4
[,1] [,2]
[1,] 1.3 0.3
[2,] 0.3 1.3

, , 5
[,1] [,2]
[1,] 1 0
[2,] 0 1
> fok(20,9,2,2)
> uu
[,1]
[1,] 201063.801383
[2,] 1302.940425
[3,] 751.572541
[4,] 230.956874
[5,] 5.566436
> ab

, , 1
[,1] [,2]
[1,] 0.4 -0.6
[2,] -0.6 0.4

, , 2
[,1] [,2]
[1,] 1.6 0.6
[2,] 0.6 1.6

, , 3
[,1] [,2]
[1,] 1.3 0.3
[2,] 0.3 1.3

, , 4
[,1] [,2]
[1,] 1 0
[2,] 0 1

, , 5
[,1] [,2]
[1,] 0.7 -0.3
[2,] -0.3 0.7
> fok(20,9,2,3)
> uu
[,1]
[1,] 1718.825011
[2,] 1302.940425
[3,] 751.572541
[4,] 230.956874
[5,] 5.566436
> ab

, , 1
[,1] [,2]
[1,] 1.9 0.9
[2,] 0.9 1.9

, , 2
[,1] [,2]
[1,] 1.6 0.6
[2,] 0.6 1.6

, , 3
[,1] [,2]
[1,] 1.3 0.3
[2,] 0.3 1.3

, , 4
[,1] [,2]
[1,] 1 0
[2,] 0 1

, , 5
[,1] [,2]
[1,] 0.7 -0.3
[2,] -0.3 0.7
> ----------------------------------------------------------------------
> fok
function(L, M, N, S)
{
fun1(L, M, N)
while(S > 0) {
S <- S - 1
fun2()
fun3(L, M, N)
}
}
> fun1
function(L, M, N)
{
ab <<- array(c(1, 0, 0, 1, 1.3, 0.3, 0.3, 1.3, 1.6, 0.6, 0.6, 1.6, 1.9,
0.9, 0.9, 1.9, 2.2, 1.2,
1.2, 2.2), dim = c(2, 2, 5))
uu <<- matrix(rep(0, 5), ncol = 1)
for(i in 1:5) {
tt <<- ab[, , i]
f2(L, M, N)
uu[i, 1] <<- u
}
}
> fun2
function()
{
ord <- rev(order(uu[, 1]))
newuu <- uu[ord, 1, drop = F]
newab <- ab[, , ord]
uu[, 1] <<- newuu
ab[, , 1:5] <<- newab
}
> fun3
function(L, M, N)
{
ab[1, 1, 1] <<- (ab[1, 1, 2] + ab[1, 1, 3] + ab[1, 1, 4] + ab[1, 1,
5])/2 - ab[1, 1, 1]
ab[1, 2, 1] <<- (ab[1, 2, 2] + ab[1, 2, 3] + ab[1, 2, 4] + ab[1, 2,
5])/2 - ab[1, 2, 1]
ab[2, 1, 1] <<- (ab[2, 1, 2] + ab[2, 1, 3] + ab[2, 1, 4] + ab[2, 1,
5])/2 - ab[2, 1, 1]
ab[2, 2, 1] <<- (ab[2, 2, 2] + ab[2, 2, 3] + ab[2, 2, 4] + ab[2, 2,
5])/2 - ab[2, 2, 1] #browser()
for(i in 1:1) {
tt <<- ab[, , i]
f2(L, M, N)
uu[i, 1] <<- u
}
#tt
}
> fun
function(L, M, N)
{
LL <- L * L
NN <- N * N
LN <- L * N
MN <- M * N
L1 <- L - 1
M1 <- M - 1
N1 <- N - 1 #---------------------------------------------------
x <- array(NA, dim = c(L, M))
for(i in 1:M)
x[, i] <- gaus[, i]
r <- matrix(rep(0, LL), ncol = L)
for(i in 1:M)
r <- r + x[, i] %*% t(x[, i])
v <- matrix(rep(0, LN), ncol = N)
v[, 1:N] <- eigen(r/9)$vectors[, 1:N]
c <- array(NA, dim = c(N, M))
c <- matrix(rep(0, MN), ncol = M)
for(i in 1:N) {
for(j in 1:M) {
for(k in 1:N) {
c[i, j] <- solve(tt)[i, k] * t(v[, k]) %*% x[,
j] + c[i, j]
}
}
}
for(i in 1:N) {
for(j in 1:M) {
if(c[i, j] > 0)
c[i, j] <- 0
else c[i, j] <- c[i, j] * c[i, j]
}
}
for(i in 1:N) {
for(j in 1:M1) {
c[i, j + 1] <- c[i, j + 1] + c[i, j]
}
}
for(i in 1:N1) {
c[i + 1, M] <- c[i + 1, M] + c[i, M]
}
penalty <- c[N, M] #----------------------------------------
s <- matrix(rep(0, LN), ncol = N)
for(i in 1:N) {
for(k in 1:N) {
s[, i] <- s[, i] + tt[i, k] * v[, k]
}
}
ds <- matrix(rep(0, LN), ncol = N)
for(j in 1:N) {
ds[1, j] <- s[2, j] - s[1, j]
ds[L, j] <- s[L, j] - s[L1, j]
for(i in 2:L1) {
ds[i, j] <- (s[i + 1, j] - s[i - 1, j])/2
}
}
dds <- matrix(rep(0, LN), ncol = N)
for(j in 1:N) {
dds[1, j] <- ds[2, j] - ds[1, j]
dds[L, j] <- ds[L, j] - ds[L1, j]
for(i in 2:L1) {
dds[i, j] <- (ds[i + 1, j] - ds[i - 1, j])/2
}
}
#----------------------------------------
abc <- matrix(rep(0, N), ncol = 1)
for(i in 1:N) {
for(k in 1:L) {
abc[i] <- abs(dds[k, i]) + abc[i]
}
}
#----------------------------------------
p <- matrix(rep(0, LN), ncol = N)
for(i in 1:N) {
for(l in 1:L) {
p[l, i] <- abs(dds[l, i])/abc[i, 1]
}
}
#----------------------------------------
h <- matrix(rep(0, N), ncol = 1)
for(i in 1:N) {
for(l in 1:L) {
h[i, 1] <- p[l, i] * log(p[l, i]) + h[i, 1]
}
}
for(i in 1:N1) {
h[i + 1, 1] <- h[i + 1, 1] + h[i, 1]
}
hi <- h[N, 1]
u <<- - hi + penalty
}
>>>>>>>>>>>>>>>>>>>>>>
"gaus"
14.06 12.54 11.02 9.50 7.98 6.47 4.95 3.43 1.91
15.72 14.06 12.40 10.73 9.07 7.41 5.75 4.09 2.43
17.05 15.31 13.57 11.84 10.10 8.36 6.62 4.88 3.15
17.96 16.23 14.49 12.76 11.03 9.30 7.57 5.83 4.10
18.37 16.74 15.10 13.47 11.84 10.21 8.58 6.94 5.31
18.25 16.81 15.36 13.92 12.48 11.04 9.60 8.16 6.71
17.61 16.43 15.24 14.06 12.88 11.70 10.52 9.34 8.16
16.49 15.61 14.72 13.83 12.95 12.06 11.18 10.29 9.41
14.98 14.38 13.79 13.19 12.59 12.00 11.40 10.80 10.20
13.18 12.83 12.47 12.12 11.77 11.41 11.06 10.71 10.35
11.22 11.04 10.86 10.68 10.50 10.32 10.14 9.97 9.79
9.22 9.14 9.07 8.99 8.91 8.83 8.76 8.68 8.60
7.32 7.28 7.24 7.21 7.17 7.13 7.09 7.05 7.01
5.60 5.57 5.53 5.49 5.46 5.42 5.38 5.35 5.31
4.14 4.09 4.04 3.99 3.94 3.88 3.83 3.78 3.73
2.95 2.89 2.82 2.76 2.69 2.63 2.56 2.50 2.43
2.04 1.97 1.90 1.83 1.76 1.69 1.62 1.55 1.48
1.36 1.30 1.23 1.17 1.10 1.03 0.97 0.90 0.84
0.88 0.83 0.77 0.72 0.66 0.61 0.55 0.50 0.45
0.55 0.51 0.47 0.43 0.39 0.35 0.31 0.27 0.22

------------------------------
Mikio
mkaihara@ichinoseki.ac.jp
Tel:+81-191-24-4772、
Fax:+81-191-24-2146
INCT, 021-8511, Japan
------------------------------
-----------------------------------------------------------------------
This message was distributed by s-news@wubios.wustl.edu. To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message: unsubscribe s-news