s-news
[Top] [All Lists]

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

To: s-news@wubios.wustl.edu
Subject: [S] Q on an optimization program (Simplex Method)
From: Mikio Kaihara <mkaihara@ichinoseki.ac.jp>
Date: Mon, 29 Jun 1998 08:22:33 +0900
Sender: owner-s-news@wubios.wustl.edu
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

<Prev in Thread] Current Thread [Next in Thread>