00001 #include "SundanceCombinatorialUtils.hpp"
00002 #include "SundanceIntVec.hpp"
00003 #include "PlayaExceptions.hpp"
00004 #include "PlayaTabs.hpp"
00005 #include "SundanceOut.hpp"
00006 #include <algorithm>
00007 #include <iterator>
00008 #include <iostream>
00009
00010
00011 namespace Sundance
00012 {
00013 using namespace Teuchos;
00014 using std::endl;
00015
00016 Array<Array<int> > partitionInteger(int n)
00017 {
00018 typedef Array<int> Aint;
00019 typedef Array<Array<int> > AAint;
00020 static Array<Array<Array<int> > > rtn =
00021 tuple<AAint>(
00022 tuple<Aint>(tuple(1)),
00023 tuple<Aint>(tuple(2), tuple(1, 1)),
00024 tuple<Aint>(
00025 tuple(3),
00026 tuple(2, 1),
00027 tuple(1, 1, 1)),
00028 tuple<Aint>(
00029 tuple(4),
00030 tuple(3, 1),
00031 tuple(2, 2),
00032 tuple(2, 1, 1),
00033 tuple(1,1,1,1)),
00034 tuple<Aint>(
00035 tuple(5),
00036 tuple(4, 1),
00037 tuple(3, 2),
00038 tuple(3, 1, 1),
00039 tuple(2, 2, 1),
00040 tuple(2, 1, 1, 1),
00041 tuple(1, 1, 1, 1, 1)
00042 ));
00043 TEUCHOS_TEST_FOR_EXCEPTION(n<1 || n>5, std::runtime_error,
00044 "case n=" << n << " not implemented in partitionInteger()");
00045 return rtn[n-1];
00046 }
00047
00048
00049 bool nextNum(Array<int>& digits, const Array<int>& radix)
00050 {
00051
00052 int n = digits.size();
00053 int j = n-1;
00054 while(j >= 0)
00055 {
00056 if (digits[j]==(radix[j]-1))
00057 {
00058 digits[j] = 0;
00059 j--;
00060 }
00061 else
00062 {
00063 digits[j]++;
00064 return true;
00065 }
00066 }
00067 return false;
00068 }
00069
00070 void weightedPartitions(int n, const Array<int>& w,
00071 Array<Array<int> >& parts)
00072 {
00073 int S = w.size();
00074 for (int i=0; i<S; i++)
00075 {
00076 TEUCHOS_TEST_FOR_EXCEPT(w[i] <= 0);
00077 }
00078
00079 Array<Array<int> > trial(S);
00080 Array<int> radix(S);
00081 for (int i=0; i<S; i++)
00082 {
00083 trial[i].resize(n/w[i]+1);
00084 for (int j=0; j<trial[i].size(); j++)
00085 {
00086 trial[i][j] = j;
00087 }
00088 radix[i] = trial[i].size();
00089 }
00090
00091 bool workLeft = true;
00092 Array<int> index(S, 0);
00093 while (workLeft)
00094 {
00095 Array<int> vals(S);
00096 int count = 0;
00097 for (int i=0; i<S; i++)
00098 {
00099 vals[i] = trial[i][index[i]];
00100 count += w[i]*vals[i];
00101 }
00102 if (count==n)
00103 {
00104 parts.append(vals);
00105 }
00106 workLeft = nextNum(index, radix);
00107 }
00108 }
00109
00110 void weightedOrderedPartitions(const IntVec& v, const Array<int>& w,
00111 Array<Array<IntVec> >& parts)
00112 {
00113 int D = v.size();
00114 int S = w.size();
00115 Array<int> radix(D);
00116
00117 Array<Array<Array<int> > > nParts(D);
00118 for (int i=0; i<D; i++)
00119 {
00120 weightedPartitions(v[i], w, nParts[i]);
00121 radix[i] = nParts[i].size();
00122 if (radix[i]==0) radix[i] = 1;
00123 }
00124
00125 bool workLeft = true;
00126 Array<int> index(D, 0);
00127 while (workLeft)
00128 {
00129 Array<IntVec> L(S, D);
00130 for (int i=0; i<D; i++)
00131 {
00132 for (int j=0; j<S; j++)
00133 {
00134 if (nParts[i].size()==0) L[i][j] = 0;
00135 else L[j][i] = nParts[i][index[i]][j];
00136 }
00137 }
00138
00139 bool accept = true;
00140 for (int i=0; i<S; i++)
00141 {
00142 if (L[i].abs()==0) {accept = false; break;}
00143 if (i>0 && !(L[i-1] < L[i])) {accept = false; break;}
00144 }
00145 if (accept)
00146 {
00147 parts.append(L);
00148 }
00149 workLeft = nextNum(index, radix);
00150 }
00151
00152 }
00153
00154
00155 void pSet(const IntVec& lambda, const IntVec& nu, int s,
00156 Array<Array<IntVec> >& K, Array<Array<IntVec> >& L)
00157 {
00158 Array<Array<IntVec> > KParts;
00159 lambda.getPartitions(s, KParts);
00160
00161 for (int i=0; i<KParts.size(); i++)
00162 {
00163 const Array<IntVec>& kPart = KParts[i];
00164 Array<int> w(s);
00165 for (int j=0; j<s; j++) w[j] = kPart[j].abs();
00166 Array<Array<IntVec> > LParts;
00167 weightedOrderedPartitions(nu, w, LParts);
00168 for (int j=0; j<LParts.size(); j++)
00169 {
00170 K.append(kPart);
00171 L.append(LParts[j]);
00172 }
00173 }
00174
00175 }
00176
00177
00178
00179 Array<Array<Array<int> > > compositions(int n)
00180 {
00181 Array<Array<Array<int> > > q(n);
00182
00183 Array<Array<int> > x = partitionInteger(n);
00184
00185 Array<Array<int> > p;
00186 for (int m=0; m<x.size(); m++)
00187 {
00188 Array<int> tmp;
00189 Array<int> y = x[m];
00190 std::sort(y.begin(), y.end());
00191 tmp.resize(y.size());
00192 std::copy(y.begin(), y.end(), tmp.begin());
00193 p.append(tmp);
00194 while (std::next_permutation(y.begin(), y.end()))
00195 {
00196 tmp.resize(y.size());
00197 std::copy(y.begin(), y.end(), tmp.begin());
00198 p.append(tmp);
00199 }
00200 }
00201
00202 for (int i=0; i<p.size(); i++)
00203 {
00204 q[p[i].size()-1].append(p[i]);
00205 }
00206
00207 return q;
00208 }
00209
00210 void restrictedCompositions(int n, int len, Array<Array<int> >& rComps)
00211 {
00212 rComps = nonNegCompositions(n, len);
00213 }
00214
00215 Array<Array<int> > nonNegCompositions(int n, int J)
00216 {
00217
00218 Array<Array<int> > parts = partitionInteger(n);
00219
00220
00221 Array<Array<int> > jParts;
00222 for (int i=0; i<parts.size(); i++)
00223 {
00224 if (parts[i].size() <= J) jParts.append(parts[i]);
00225 }
00226
00227
00228 for (int i=0; i<jParts.size(); i++)
00229 {
00230 for (int j=jParts[i].size(); j<J; j++)
00231 {
00232 jParts[i].append(0);
00233 }
00234 }
00235
00236
00237 Array<Array<int> > all;
00238 for (int i=0; i<jParts.size(); i++)
00239 {
00240 Array<int> tmp;
00241 Array<int> y = jParts[i];
00242 std::sort(y.begin(), y.end());
00243 tmp.resize(y.size());
00244 std::copy(y.begin(), y.end(), tmp.begin());
00245 all.append(tmp);
00246 while (std::next_permutation(y.begin(), y.end()))
00247 {
00248 tmp.resize(y.size());
00249 std::copy(y.begin(), y.end(), tmp.begin());
00250 all.append(tmp);
00251 }
00252 }
00253
00254 return all;
00255 }
00256
00257
00258
00259 Set<Pair<MultiSet<int> > >
00260 loadPartitions(int x, int n,
00261 const MultiSet<int>& left,
00262 const MultiSet<int>& right)
00263 {
00264 Set<Pair<MultiSet<int> > > rtn;
00265 for (int i=0; i<n; i++)
00266 {
00267 MultiSet<int> L = left;
00268 MultiSet<int> R = right;
00269
00270 for (int j=0; j<n; j++)
00271 {
00272 if (j < i) L.put(x);
00273 else R.put(x);
00274 }
00275 rtn.put(Pair<MultiSet<int> >(L, R));
00276 rtn.put(Pair<MultiSet<int> >(R, L));
00277 }
00278 return rtn;
00279 }
00280
00281
00282 Set<Pair<MultiSet<int> > >
00283 binaryPartition(const MultiSet<int>& m)
00284 {
00285 Set<Pair<MultiSet<int> > > tmp1;
00286
00287
00288 Map<int, int> counts = countMap(m);
00289
00290 for (Map<int, int>::const_iterator
00291 i=counts.begin(); i!=counts.end(); i++)
00292 {
00293 int x = i->first;
00294 int n = i->second;
00295 if (tmp1.size()==0)
00296 {
00297 MultiSet<int> L;
00298 MultiSet<int> R;
00299 tmp1 = loadPartitions(x, n, L, R);
00300 }
00301 else
00302 {
00303 Set<Pair<MultiSet<int> > > tmp2;
00304 for (Set<Pair<MultiSet<int> > >::const_iterator
00305 j=tmp1.begin(); j!=tmp1.end(); j++)
00306 {
00307 MultiSet<int> L = j->first();
00308 MultiSet<int> R = j->second();
00309 Set<Pair<MultiSet<int> > > t
00310 = loadPartitions(x, n, L, R);
00311 tmp2.merge(t);
00312 }
00313 tmp1 = tmp2;
00314 }
00315 }
00316
00317 #ifdef BLAH
00318 Set<SortedPair<MultiSet<int> > > rtn;
00319 for (Set<Pair<MultiSet<int> > >::const_iterator
00320 j=tmp1.begin(); j!=tmp1.end(); j++)
00321 {
00322 rtn.put(SortedPair<MultiSet<int> >(j->first(), j->second()));
00323 }
00324 #endif
00325
00326
00327 return tmp1;
00328 }
00329
00330
00331
00332 Set<MultiSet<MultiSet<int> > >
00333 multisetPartitions(const MultiSet<int>& m)
00334 {
00335 Set<MultiSet<MultiSet<int> > > rtn;
00336 int mSize = m.size();
00337 if (m.size()==0) return rtn;
00338 if (m.size()==1) return makeSet(makeMultiSet(m));
00339 rtn.put(makeMultiSet(m));
00340
00341 Set<Pair<MultiSet<int> > >
00342 twoParts = binaryPartition(m);
00343
00344 for (Set<Pair<MultiSet<int> > >::const_iterator
00345 i=twoParts.begin(); i!=twoParts.end(); i++)
00346 {
00347 MultiSet<int> L = i->first();
00348 MultiSet<int> R = i->second();
00349 if ((int) L.size() == mSize || (int) R.size() == mSize) continue;
00350 if ((int) L.size() > 0 && (int) R.size() > 0) rtn.put(makeMultiSet(L, R));
00351
00352 Array<MultiSet<MultiSet<int> > > lParts;
00353 Array<MultiSet<MultiSet<int> > > rParts;
00354 if (L.size() > 0)
00355 {
00356 lParts = multisetPartitions(L).elements();
00357 }
00358 if (R.size() > 0)
00359 {
00360 rParts = multisetPartitions(R).elements();
00361 }
00362 for (int j=0; j<lParts.size(); j++)
00363 {
00364 for (int k=0; k<rParts.size(); k++)
00365 {
00366 const MultiSet<MultiSet<int> >& a = lParts[j];
00367 const MultiSet<MultiSet<int> >& b = rParts[k];
00368 MultiSet<MultiSet<int> > combined;
00369 for (MultiSet<MultiSet<int> >::const_iterator
00370 x=a.begin(); x!=a.end(); x++)
00371 {
00372 combined.put(*x);
00373 }
00374 for (MultiSet<MultiSet<int> >::const_iterator
00375 x=b.begin(); x!=b.end(); x++)
00376 {
00377 combined.put(*x);
00378 }
00379 rtn.put(combined);
00380 }
00381 }
00382 }
00383 return rtn;
00384 }
00385
00386
00387 Map<int, int> countMap(const MultiSet<int>& m)
00388 {
00389 Map<int, int> rtn;
00390 for (MultiSet<int>::const_iterator i=m.begin(); i!=m.end(); i++)
00391 {
00392 if (rtn.containsKey(*i)) rtn[*i]++;
00393 else rtn.put(*i, 1);
00394 }
00395 return rtn;
00396 }
00397
00398
00399 Array<Array<int> > indexCombinations(const Array<int>& s)
00400 {
00401
00402 Array<int> D(s.size(), 0);
00403 Array<int> C(s.size(), -1);
00404
00405 int p=1;
00406 for (int k=1; k<=s.size(); k++)
00407 {
00408 D[s.size()-k] = p;
00409 p = p*s[s.size() - k];
00410 }
00411 Array<Array<int> > rtn(p);
00412
00413 for (int i=0; i<p; i++)
00414 {
00415 for (int j=0; j<s.size(); j++)
00416 {
00417 if ( (i % D[j])==0 )
00418 {
00419 C[j] = (C[j]+1) % s[j];
00420 }
00421 }
00422 rtn[i] = C;
00423 }
00424 return rtn;
00425 }
00426
00427
00428 Array<Array<Array<int> > > binnings(const MultiSet<int>& mu, int n)
00429 {
00430 int N = mu.size();
00431 Array<Array<int> > c = compositions(N)[n-1];
00432 Array<Array<Array<int> > > rtn;
00433
00434 for (int i=0; i<c.size(); i++)
00435 {
00436 Array<Array<Array<int> > > a = indexArrangements(mu, c[i]);
00437 for (int j=0; j<a.size(); j++)
00438 {
00439 rtn.append(a[j]);
00440 }
00441 }
00442 return rtn;
00443 }
00444
00445
00446
00447 Array<Array<Array<int> > > indexArrangements(const MultiSet<int>& mu,
00448 const Array<int> & k)
00449 {
00450 int nBins = k.size();
00451
00452 int M = 0;
00453 for (int i=0; i<nBins; i++)
00454 {
00455 M += k[i];
00456 }
00457
00458 Array<int> I;
00459 for (MultiSet<int>::const_iterator iter=mu.begin(); iter!=mu.end(); iter++)
00460 {
00461 I.append(*iter);
00462 }
00463
00464 Array<Array<Array<int> > > rtn;
00465
00466 do
00467 {
00468 Array<Array<int> > bins(nBins);
00469 int count = 0;
00470 for (int i=0; i<nBins; i++)
00471 {
00472 for (int j=0; j<k[i]; j++)
00473 {
00474 bins[i].append(I[count++]);
00475 }
00476 }
00477 rtn.append(bins);
00478 }
00479 while (std::next_permutation(I.begin(), I.end()));
00480 return rtn;
00481
00482 }
00483
00484
00485 Set<MultiSet<int> > multisetSubsets(const MultiSet<int>& x)
00486 {
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498 Array<int> elements = x.elements();
00499
00500
00501
00502 int n = elements.size();
00503 int maxNumSubsets = pow2(n);
00504
00505 Set<MultiSet<int> > rtn;
00506
00507
00508
00509
00510 for (int i=1; i<maxNumSubsets; i++)
00511 {
00512 Array<int> bits = bitsOfAnInteger(i, n);
00513 MultiSet<int> ms;
00514 for (int j=0; j<n; j++)
00515 {
00516 if (bits[j] == 1) ms.put(elements[j]);
00517 }
00518 rtn.put(ms);
00519 }
00520 return rtn;
00521 }
00522
00523
00524 Array<Array<MultiSet<int> > > multisetCompositions(int s,
00525 const MultiSet<int>& x)
00526 {
00527 Set<MultiSet<MultiSet<int> > > parts = multisetPartitions(x);
00528
00529 Array<Array<MultiSet<int> > > rtn;
00530
00531 typedef Set<MultiSet<MultiSet<int> > >::const_iterator iter;
00532 for (iter i=parts.begin(); i!=parts.end(); i++)
00533 {
00534 if ((int) i->size()!=s) continue;
00535 Array<MultiSet<int> > y = i->elements();
00536 Array<MultiSet<int> > tmp(y.size());
00537 std::sort(y.begin(), y.end());
00538 std::copy(y.begin(), y.end(), tmp.begin());
00539 rtn.append(tmp);
00540 while (std::next_permutation(y.begin(), y.end()))
00541 {
00542 tmp.resize(y.size());
00543 std::copy(y.begin(), y.end(), tmp.begin());
00544 rtn.append(tmp);
00545 }
00546 }
00547 return rtn;
00548 }
00549
00550 Array<Array<int> > distinctIndexTuples(int m, int n)
00551 {
00552 Array<Array<int> > rtn;
00553 Array<int> cur(m);
00554 for (int i=0; i<m; i++) cur[i] = i;
00555 if (m==1)
00556 {
00557 for (int i=0; i<n; i++) rtn.append(tuple(i));
00558 return rtn;
00559 }
00560
00561 while(cur[0] < (n-m)+1)
00562 {
00563 rtn.append(cur);
00564 for (int i=(m-1); i>=0; i--)
00565 {
00566 if (cur[i] < (n-1))
00567 {
00568 cur[i]++;
00569 bool overflow = false;
00570 for (int j=i+1; j<m; j++)
00571 {
00572 cur[j] = cur[j-1]+1;
00573 if (cur[j] >= n) overflow = true;
00574 }
00575 if (overflow) continue;
00576 else break;
00577 }
00578 }
00579 }
00580 return rtn;
00581 }
00582
00583 Array<int> bitsOfAnInteger(int x, int n)
00584 {
00585 TEUCHOS_TEST_FOR_EXCEPTION(x >= pow2(n), std::logic_error,
00586 "Invalid input to bitsOfAnIteger");
00587
00588 Array<int> rtn(n);
00589
00590 int r = x;
00591 for (int b=n-1; b>=0; b--)
00592 {
00593 rtn[b] = r/pow2(b);
00594 r = r - rtn[b]*pow2(b);
00595 }
00596 return rtn;
00597 }
00598
00599
00600
00601
00602
00603 int pow2(int n)
00604 {
00605 static Array<int> p2(1,1);
00606
00607 if ((int) n >= p2.size())
00608 {
00609 int oldN = p2.size();
00610 for (int i=oldN; i<=n; i++) p2.push_back(p2[i-1]*2);
00611 }
00612
00613 return p2[n];
00614 }
00615 }
00616
00617
00618
00619
00620
00621
00622
00623
00624