optimization challenge

FYI, in bioinformatics, we use dynamic programming algorithms in similar
ways to solve similar problems of finding guaranteed-optimal partitions in
streams of data (usually DNA or protein sequence, but sometimes numerical
data from chip-arrays).  These "path optimization" algorithms are often
called Viterbi algorithms, a web search for which should provide multiple
references.

The solutions are not necessarily unique (there may be multiple
paths/partitions with identical integer maxima in some systems) and there is
much research on whether the optimal solution is actually the one you want
to work with (for example, there may be a fair amount of probability mass
within an area/ensemble of suboptimal solutions that overall have greater
posterior probabilities than does the optimal solution "singleton").  See
Chip Lawrence's PNAS paper for more erudite discussion, and references
therein: www.pnas.org/content/105/9/3209.abstract

-Aaron

P.S. Good to see you here Albyn -- I enjoyed your stat. methods course at
Reed back in 1993, which started me down a somewhat windy road to
statistical genomics!

--
Aaron J. Mackey, PhD
Assistant Professor
Center for Public Health Genomics
University of Virginia
amackey at virginia.edu


On Wed, Jan 13, 2010 at 5:23 PM, Ravi Varadhan <RVaradhan at jhmi.edu> wrote:

Greg - thanks for posting this interesting problem.

Albyn - thanks for posting a solution.  Now, I have some questions: (1) is
the algorithm guaranteed to find a "best" solution? (2) can there be
multiple solutions (it seems like there can be more than 1 solution
depending on the data)?, and (3) is there a good reference for this and
similar algorithms?

Thanks & Best,
Ravi.


----------------------------------------------------------------------------
-------

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: rvaradhan at jhmi.edu

Webpage:

http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h
tml<http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h%0Atml>




----------------------------------------------------------------------------
--------


-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
On
Behalf Of Albyn Jones
Sent: Wednesday, January 13, 2010 1:19 PM
To: Greg Snow
Cc: r-help at r-project.org
Subject: Re: [R] optimization challenge

The key idea is that you are building a matrix that contains the
solutions to smaller problems which are sub-problems of the big
problem.  The first row of the matrix SSQ contains the solution for no
splits, ie SSQ[1,j] is just the sum of squares about the overall mean
for reading chapters1 through j in one day.  The iteration then uses
row m-1 to construct row m, since if SSQ[m-1,j] (optimal reading of j
chapters in m-1 days) is part of the overall optimal solution, you
have already computed it, and so don't ever need to recompute it.

      TS = SSQ[m-1,j]+(SSQ1[j+1])

computes the vector of possible solutions for SSQ[m,n] (n chapters in n
days)
breaking it into two pieces: chapters 1 to j in m-1 days, and chapters j+1
to
n in 1 day.  j is a vector in the function, and min(TS) is the minimum
over choices of j, ie SSQ[m,n].

At the end, SSQ[128,239] is the optimal value for reading all 239
chapters in 128 days.  That's just the objective function, so the rest
involves constructing the list of optimal cuts, ie which chapters are
grouped together for each day's reading.  That code uses the same
idea... constructing a list of lists of cutpoints.

statisticians should study a bit of data structures and algorithms!

albyn

On Wed, Jan 13, 2010 at 10:45:11AM -0700, Greg Snow wrote:

WOW, your results give about half the variance of my best optim run

(possibly due to my suboptimal use of optim).

Can you describe a little what the algorithm is doing?

--
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.snow at imail.org
801.408.8111

-----Original Message-----
From: Albyn Jones [mailto:jones at reed.edu]
Sent: Tuesday, January 12, 2010 5:31 PM
To: Greg Snow
Cc: r-help at r-project.org
Subject: Re: [R] optimization challenge

Greg

Nice problem: I wasted my whole day on it :-)

I was explaining my plan for a solution to a colleague who is a
computer scientist, he pointed out that I was trying to re-invent the
wheel known as dynamic programming.  here is my code, apparently it is
called "bottom up dynamic programming".  It runs pretty quickly, and
returns (what I hope is :-) the optimal sum of squares and the
cut-points.

function(X=bom3$Verses,days=128){
# find optimal BOM reading schedule for Greg Snow
# minimize variance of quantity to read per day over 128 days
#
N = length(X)
Nm1 = N-1
SSQ<- matrix(NA,nrow=days,ncol=N)
Cuts <- list()
#
#  SSQ[i,j]: the ssqs about the overall mean for the optimal partition
#           for i days on the chapters 1 to j
#
M = sum(X)/days
CS = cumsum(X)
SSQ[1,]= (CS-M)^2
Cuts[[1]]= as.list(1:N)
#
for(m in 2:days){
Cuts[[m]]=list()
#for(i in 1:(m-1)) Cuts[[m]][[i]] = Cuts[[m-1]][[i]]
for(n in m:N){
                      CS = cumsum(X[n:1])[n:1]
                      SSQ1 = (CS-M)^2
              j = (m-1):(n-1)
                      TS = SSQ[m-1,j]+(SSQ1[j+1])
              SSQ[m,n] = min(TS)
                  k = min(which((min(TS)== TS)))+m-1
                  Cuts[[m]][[n]] = c(Cuts[[m-1]][[k-1]],n)
}
}
list(SSQ=SSQ[days,N],Cuts=Cuts[[days]][[N]])
}

$SSQ
[1] 11241.05

$Cuts
  [1]   2   4   7   9  11  13  15  16  17  19  21  23  25  27  30  31
34  37
 [19]  39  41  44  46  48  50  53  56  59  60  62  64  66  68  70  73
75  77
 [37]  78  80  82  84  86  88  89  91  92  94  95  96  97  99 100 103
105 106
 [55] 108 110 112 113 115 117 119 121 124 125 126 127 129 131 132 135
137 138
 [73] 140 141 142 144 145 146 148 150 151 152 154 156 157 160 162 163
164 166
 [91] 167 169 171 173 175 177 179 181 183 185 186 188 190 192 193 194
196 199
[109] 201 204 205 207 209 211 213 214 215 217 220 222 223 225 226 228
234 236
[127] 238 239




On Tue, Jan 12, 2010 at 11:33:36AM -0700, Greg Snow wrote:

I have a challenge that I want to share with the group.

This is not homework (but I may assign it as such if I teach the

appropriate class again) and I have found one solution, so don't need
anything urgent.  This is more for fun to see if others can find a
better solution than I did.

The challenge:

I want to read a book in a given number of days.  I want to read an

integer number of chapters each day (there are more chapters than
days), no stopping part way through a chapter, and at least 1 chapter
each day.  The chapters are very non uniform in length (some very
short, a few very long, many in between) so I would like to come up
with a reading schedule that minimizes the variance of the length of
the days readings (read multiple short chapters on the same day, long
chapters are the only one read that day).  I also want to read through
the book in order (no skipping ahead to combine short chapters that are
not naturally next to each other.

My thought was that the optim function with method="SANN" would be an

appropriate approach, but my first couple of tries did not give very
good results.  I have since come up with an optim with SANN solution
that gives what I consider good results (but I accept that better is
possible).

Below is a data frame with the lengths of the chapters for the book

that originally sparked the challenge for me (but the general idea
should work for any book).  Each row represents a chapter (in order)
with 3 different measures of the length of the chapter.

For this challenge I want to read the book in 128 days (there are 239

chapters).

I will post my solutions in a few days, but I want to wait so that my

direction does not influence people from trying other approaches (if
there is something better than optim, that is fine).

Good luck for anyone interested in the challenge,

The data frame:

bom3 <- structure(list(Chapter = structure(1:239, .Label = c("1 Nephi

1",

"1 Nephi 2", "1 Nephi 3", "1 Nephi 4", "1 Nephi 5", "1 Nephi 6",
"1 Nephi 7", "1 Nephi 8", "1 Nephi 9", "1 Nephi 10", "1 Nephi 11",
"1 Nephi 12", "1 Nephi 13", "1 Nephi 14", "1 Nephi 15", "1 Nephi 16",
"1 Nephi 17", "1 Nephi 18", "1 Nephi 19", "1 Nephi 20", "1 Nephi 21",
"1 Nephi 22", "2 Nephi 1", "2 Nephi 2", "2 Nephi 3", "2 Nephi 4",
"2 Nephi 5", "2 Nephi 6", "2 Nephi 7", "2 Nephi 8", "2 Nephi 9",
"2 Nephi 10", "2 Nephi 11", "2 Nephi 12", "2 Nephi 13", "2 Nephi 14",
"2 Nephi 15", "2 Nephi 16", "2 Nephi 17", "2 Nephi 18", "2 Nephi 19",
"2 Nephi 20", "2 Nephi 21", "2 Nephi 22", "2 Nephi 23", "2 Nephi 24",
"2 Nephi 25", "2 Nephi 26", "2 Nephi 27", "2 Nephi 28", "2 Nephi 29",
"2 Nephi 30", "2 Nephi 31", "2 Nephi 32", "2 Nephi 33", "Jacob 1",
"Jacob 2", "Jacob 3", "Jacob 4", "Jacob 5", "Jacob 6", "Jacob 7",
"Enos 1", "Jarom 1", "Omni 1", "Words of Mormon 1", "Mosiah 1",
"Mosiah 2", "Mosiah 3", "Mosiah 4", "Mosiah 5", "Mosiah 6", "Mosiah

7",

"Mosiah 8", "Mosiah 9", "Mosiah 10", "Mosiah 11", "Mosiah 12",
"Mosiah 13", "Mosiah 14", "Mosiah 15", "Mosiah 16", "Mosiah 17",
"Mosiah 18", "Mosiah 19", "Mosiah 20", "Mosiah 21", "Mosiah 22",
"Mosiah 23", "Mosiah 24", "Mosiah 25", "Mosiah 26", "Mosiah 27",
"Mosiah 28", "Mosiah 29", "Alma 1", "Alma 2", "Alma 3", "Alma 4",
"Alma 5", "Alma 6", "Alma 7", "Alma 8", "Alma 9", "Alma 10",
"Alma 11", "Alma 12", "Alma 13", "Alma 14", "Alma 15", "Alma 16",
"Alma 17", "Alma 18", "Alma 19", "Alma 20", "Alma 21", "Alma 22",
"Alma 23", "Alma 24", "Alma 25", "Alma 26", "Alma 27", "Alma 28",
"Alma 29", "Alma 30", "Alma 31", "Alma 32", "Alma 33", "Alma 34",
"Alma 35", "Alma 36", "Alma 37", "Alma 38", "Alma 39", "Alma 40",
"Alma 41", "Alma 42", "Alma 43", "Alma 44", "Alma 45", "Alma 46",
"Alma 47", "Alma 48", "Alma 49", "Alma 50", "Alma 51", "Alma 52",
"Alma 53", "Alma 54", "Alma 55", "Alma 56", "Alma 57", "Alma 58",
"Alma 59", "Alma 60", "Alma 61", "Alma 62", "Alma 63", "Helaman 1",
"Helaman 2", "Helaman 3", "Helaman 4", "Helaman 5", "Helaman 6",
"Helaman 7", "Helaman 8", "Helaman 9", "Helaman 10", "Helaman 11",
"Helaman 12", "Helaman 13", "Helaman 14", "Helaman 15", "Helaman 16",
"3 Nephi 1", "3 Nephi 2", "3 Nephi 3", "3 Nephi 4", "3 Nephi 5",
"3 Nephi 6", "3 Nephi 7", "3 Nephi 8", "3 Nephi 9", "3 Nephi 10",
"3 Nephi 11", "3 Nephi 12", "3 Nephi 13", "3 Nephi 14", "3 Nephi 15",
"3 Nephi 16", "3 Nephi 17", "3 Nephi 18", "3 Nephi 19", "3 Nephi 20",
"3 Nephi 21", "3 Nephi 22", "3 Nephi 23", "3 Nephi 24", "3 Nephi 25",
"3 Nephi 26", "3 Nephi 27", "3 Nephi 28", "3 Nephi 29", "3 Nephi 30",
"4 Nephi 1", "Mormon 1", "Mormon 2", "Mormon 3", "Mormon 4",
"Mormon 5", "Mormon 6", "Mormon 7", "Mormon 8", "Mormon 9", "Ether

1",

"Ether 2", "Ether 3", "Ether 4", "Ether 5", "Ether 6", "Ether 7",
"Ether 8", "Ether 9", "Ether 10", "Ether 11", "Ether 12", "Ether 13",
"Ether 14", "Ether 15", "Moroni 1", "Moroni 2", "Moroni 3", "Moroni

4",

"Moroni 5", "Moroni 6", "Moroni 7", "Moroni 8", "Moroni 9", "Moroni

10"

), class = "factor"), Words = c(908L, 879L, 1067L, 1262L, 761L,
202L, 992L, 1221L, 259L, 924L, 1315L, 860L, 1899L, 1284L, 1488L,
1618L, 2523L, 1217L, 1292L, 698L, 945L, 1506L, 1543L, 1460L,
1170L, 1300L, 1169L, 895L, 405L, 812L, 2388L, 966L, 338L, 647L,
587L, 203L, 857L, 370L, 687L, 570L, 587L, 928L, 520L, 134L, 587L,
891L, 1699L, 1483L, 1461L, 1240L, 804L, 708L, 988L, 426L, 647L,
719L, 1365L, 619L, 929L, 3758L, 511L, 1242L, 1160L, 734L, 1398L,
857L, 966L, 2112L, 1117L, 1605L, 740L, 309L, 1555L, 938L, 864L,
957L, 1271L, 1291L, 1073L, 391L, 1056L, 560L, 650L, 1382L, 978L,
963L, 1328L, 620L, 1186L, 954L, 837L, 1200L, 1601L, 812L, 1879L,
1496L, 1433L, 1016L, 1035L, 2787L, 390L, 1443L, 1231L, 1511L,
1442L, 1466L, 1858L, 1347L, 1526L, 711L, 1010L, 1848L, 1623L,
1701L, 1279L, 955L, 1871L, 764L, 1509L, 752L, 1665L, 1260L, 575L,
708L, 2666L, 1478L, 1837L, 855L, 1576L, 699L, 1229L, 2027L, 651L,
793L, 1153L, 701L, 1229L, 2221L, 1243L, 911L, 1809L, 1572L, 1073L,
1345L, 1734L, 1682L, 1757L, 1085L, 988L, 1218L, 2177L, 1517L,
1748L, 489L, 1777L, 854L, 2220L, 593L, 1418L, 599L, 1527L, 1076L,
2084L, 1797L, 1142L, 1210L, 1480L, 720L, 1561L, 873L, 1855L,
1241L, 824L, 1025L, 1340L, 769L, 1353L, 1518L, 931L, 1200L, 1132L,
871L, 965L, 807L, 1434L, 1294L, 834L, 628L, 731L, 913L, 879L,
1344L, 1233L, 1538L, 1155L, 507L, 415L, 620L, 183L, 793L, 1269L,
1457L, 394L, 130L, 1949L, 706L, 1262L, 939L, 822L, 1067L, 914L,
454L, 1710L, 1510L, 901L, 1370L, 1253L, 892L, 226L, 1048L, 899L,
1231L, 1424L, 1415L, 762L, 1536L, 1219L, 1132L, 1314L, 147L,
119L, 120L, 153L, 91L, 335L, 1929L, 1064L, 1012L, 1149L), Letters =

c(3746L,

3640L, 4295L, 4968L, 3202L, 777L, 3941L, 4769L, 1075L, 3829L,
5159L, 3527L, 7874L, 5236L, 6149L, 6555L, 10180L, 4960L, 5399L,
2767L, 3758L, 6275L, 6391L, 6151L, 4685L, 5267L, 4767L, 3780L,
1574L, 3261L, 9985L, 4106L, 1332L, 2557L, 2454L, 810L, 3540L,
1406L, 2668L, 2304L, 2474L, 3715L, 2090L, 533L, 2442L, 3587L,
7294L, 6363L, 5809L, 4997L, 3197L, 2941L, 4065L, 1743L, 2539L,
3149L, 5815L, 2765L, 4027L, 15366L, 2091L, 4988L, 4730L, 3149L,
5906L, 3704L, 4188L, 8675L, 4771L, 6586L, 2965L, 1313L, 6326L,
3948L, 3489L, 3942L, 5205L, 5309L, 4300L, 1574L, 4559L, 2393L,
2655L, 5755L, 4070L, 4027L, 5649L, 2651L, 4978L, 3970L, 3610L,
5112L, 6733L, 3513L, 7990L, 6356L, 6179L, 4340L, 4432L, 11260L,
1657L, 5913L, 5033L, 6274L, 5835L, 5865L, 7802L, 5776L, 6295L,
3034L, 4383L, 7842L, 6659L, 6905L, 5104L, 4162L, 7733L, 3321L,
6241L, 3212L, 6840L, 5294L, 2471L, 2748L, 10762L, 6269L, 7479L,
3416L, 6477L, 3043L, 4719L, 8764L, 2494L, 3168L, 4790L, 2983L,
5129L, 9795L, 5043L, 3913L, 7635L, 6693L, 4669L, 5968L, 7508L,
7271L, 7543L, 4730L, 4095L, 5152L, 9129L, 6395L, 7357L, 2166L,
7512L, 3544L, 9554L, 2505L, 6171L, 2568L, 6542L, 4781L, 8655L,
7721L, 4751L, 5134L, 6032L, 3013L, 6536L, 3504L, 7527L, 5106L,
3607L, 4247L, 5524L, 3404L, 5947L, 6506L, 3891L, 5256L, 4944L,
3700L, 4006L, 3392L, 5732L, 5185L, 3402L, 2451L, 2962L, 3679L,
3550L, 5429L, 5083L, 6354L, 4659L, 2155L, 1754L, 2496L, 721L,
3326L, 4917L, 5856L, 1589L, 564L, 8300L, 2987L, 5317L, 3820L,
3595L, 4483L, 3664L, 1805L, 6931L, 6225L, 3508L, 5481L, 5007L,
3623L, 891L, 4265L, 3746L, 5041L, 5823L, 5726L, 3238L, 6443L,
5053L, 4664L, 5370L, 614L, 459L, 491L, 629L, 351L, 1457L, 7812L,
4677L, 4305L, 4603L), Verses = c(20, 24, 31, 38, 22, 6, 22, 38,
6, 22, 36, 23, 42, 30, 36, 39, 55, 25, 24, 22, 26, 31, 32, 30,
25, 35, 34, 18, 11, 25, 54, 25, 8, 22, 26, 6, 30, 13, 25, 22,
21, 34, 16, 6, 22, 32, 30, 33, 35, 32, 14, 18, 21, 9, 15, 19,
35, 14, 18, 77, 13, 27, 27, 15, 30, 18, 18, 41, 27, 30, 15, 7,
33, 21, 19, 22, 29, 37, 35, 12, 31, 15, 20, 35, 29, 26, 36, 16,
39, 25, 24, 39, 37, 20, 47, 33, 38, 27, 20, 62, 8, 27, 32, 34,
32, 46, 37, 31, 29, 19, 21, 39, 43, 36, 30, 23, 35, 18, 30, 17,
37, 30, 14, 17, 60, 38, 43, 23, 41, 16, 30, 47, 15, 19, 26, 15,
31, 54, 24, 24, 41, 36, 25, 30, 40, 37, 40, 23, 24, 35, 57, 36,
41, 13, 36, 21, 52, 17, 34, 14, 37, 26, 52, 41, 29, 28, 41, 19,
38, 26, 39, 31, 17, 25, 30, 19, 26, 33, 26, 30, 26, 25, 22, 19,
41, 48, 34, 27, 24, 20, 25, 39, 36, 46, 29, 17, 14, 18, 6, 21,
33, 39, 9, 2, 49, 19, 29, 22, 23, 24, 22, 10, 41, 37, 43, 25,
28, 19, 6, 30, 27, 26, 35, 34, 23, 41, 31, 31, 34, 4, 3, 4, 3,
2, 9, 48, 30, 26, 34)), .Names = c("Chapter", "Words", "Letters",
"Verses"), row.names = c(NA, -239L), class = "data.frame")


--
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.snow at imail.org
801.408.8111

guide.html