share/cobyla/ex/transport.mac

   1 /*
   2  * Maximum entropy estimate for a transportation problem.  Based on an
   3  * example from Kostas Oikonomou.
   4  *
   5  * There are 5 cities connected by two highways:
   6  *
   7  *  1 ----- 2 ----- 3 ----- 4
   8  *  |               |
   9  *  5---------------+
  10  *
  11  * So cities 1,2,3,4 are connected by one highway, and there's another
  12  * highway from 1 to 5 to 3.
  13  *
  14  * The number of cars the left cities 1, 2, and 4 is known.  The
  15  * number of cars that travelled the highway segment from 2 to 3 is
  16  * also known, and finally it is known that at least a certain number
  17  * travelled the 5 to 3 segment.  Given this information, we want to
  18  * know the most likely number (fraction) of cars that travelled
  19  * between each pair of cities on that day.  This is the 5x5 matrix C
  20  * = [c[i,j]] where c[i,i] represents the fraction of cars that left
  21  * city i and return to it.  Our information is:
  22  *
  23  * sum(c[i,j], j, 1, 5) = r[i] for i = 1, 2, 4
  24  *
  25  * c[1,3] + c[1,4] + c[2,4] = s[2,3] and
  26  * c[1,3] + c[5,3] + c[1,4] >= s[5,3].
  27  *
  28  * Suppose (r1,r2,r4) = (.13,.25,.1) and (s23,s53) = (0.11, .07).
  29  * Then the maximum entropy solution is
  30  *
  31  * matrix([.030790, .030790, .018816, .018816, .030790],
  32  *        [.059210, .059210, .059210, .059210, .059210],
  33  *        [.052,    .052,    .052,    .052,    .052   ],
  34  *        [.02,     .02,     .02,     .02,     .02    ],
  35  *        [.052,    .052,    .052,    .052,    .052   ],
  36  */
  37
  38 HH(X) := -xreduce("+",map(lambda([x],x*log(x)),X));
  39
  40 C :
  41 [c11,c12,c13,c14,c15,
  42 c21,c22,c23,c24,c25,
  43 c31,c32,c33,c34,c35,
  44 c41,c42,c43,c44,c45,
  45 c51,c52,c53,c54,c55]$
  46
  47 cars(r,s) :=
  48   fmin_cobyla(-HH(C) ,
  49               C,
  50               makelist(1/25,k,1,25),
  51               constraints =
  52                [xreduce("+",C)=1,
  53                 c11+c12+c13+c14+c15=r[1],
  54                 c21+c22+c23+c24+c25=r[2],
  55                 c41+c42+c43+c44+c45=r[3],
  56                 c13+c14+c23+c24=s[1],
  57                 c13+c53+c14>=s[2]
  58                 ],
  59               iprint=1,
  60               rhobeg = .1,
  61               rhoend = 1d-9
  62 );
  63
  64 cars([.13,.25,.1],[.11,.07]);
  65 /*
  66
  67    Normal return from subroutine COBYLA
  68
  69    NFVALS =  944   F =-3.141900E+00    MAXCV = 1.082467E-15
  70    X = 3.078947E-02   3.078947E-02   1.881578E-02   1.881579E-02   3.078948E-02
  71        5.921053E-02   5.921052E-02   3.618421E-02   3.618421E-02   5.921053E-02
  72    5.199999E-02   5.200000E-02   5.200001E-02   5.200001E-02
  73    5.200000E-02   2.000001E-02   2.000000E-02   2.000000E-02
  74    2.000000E-02   2.000000E-02   5.200000E-02   5.199999E-02
  75    5.199998E-02   5.200000E-02   5.200001E-02
  76
  77 (%o19) [[c11 = .03078947216080049, c12 = .03078947258983758,
  78 c13 = .01881578424530663, c14 = .01881579235550251, c15 = .03078947864855278,
  79 c21 = .05921052755891571, c22 = .05921051666669622, c23 = .03618421394514965,
  80 c24 = .03618420945404111, c25 = .05921053237519735, c31 = .05199998745697864,
  81 c32 = .05199999998006931, c33 = .05200001440859302, c34 = .05200000844126493,
  82 c35 = .05200000226927205, c41 = .02000000599466734, c42 = .01999999747226151,
  83 c43 = .01999999589953518, c44 = .02000000095933785, c45 = .01999999967419817,
  84 c51 = .05200000159627225, c52 = .05199999448887869, c53 = 0.0519999846606873,
  85 c54 = .05199999603531405, c55 = .05200001066267089], - 3.141900424608619,
  86 944, 0]
  87 */
  88