#include #include #include #include #include using namespace std; // Returns value of Binomial Coefficient C(n, k) int binomial(int n, int k) { int res = 1; if (n<=0 || k<=0 || k>=n) return res; if ( k > n - k ) k = n - k; for (int i = 0; i < k; ++i) { res *= (n - i); res /= (i + 1); } return res; } // Computes the fraction between perm(a') and perm(a) double prob_perm(vector flats, vector flats_prime, vector peaks, vector peaks_prime) { int h1 = flats.size(); double p_perm = 1; for (int i=1; i peaks(1,0); vector flats; vector a; int h,h_max,v,stepnumber,seed; // READ INPUT cin >> h >> h_max; cin >> stepnumber; cin >> seed; for (int k=0; k<2*h+1; k++) { cin >> v; if (k%2==0) {flats.push_back(v);} else {peaks.push_back(v);} } peaks.push_back(0); flats.push_back(0); // START RUNNING THE MARKOVIAN PROCESS vector peaks_prime; vector flats_prime; int i; int j; int op; // INIT RANDOM GENERATORS mt19937 gen(seed); uniform_real_distribution randprob(0,1); uniform_int_distribution randindex(0, h_max); uniform_int_distribution randop(1, 8); while(stepnumber>0) { stepnumber-=1; // SELECT 'i' AND 'j' IN {0,1,...,h_max} AT RANDOM. // SELECT OPERATION 'op' IN {1,2,3,4,5,6,7,8} AT RANDOM. i = randindex(gen); j = randindex(gen); op = randop(gen); peaks_prime = peaks; flats_prime = flats; if (i>h+1 || j>h+1) continue; // APPLY 'op(i,j)' TO 'a' AND OBTAIN 'a_prime' IF STEP IS FEASIBLE switch (op) { case 1: // PF(i,j) // TEST IF FEASIBLE AND APPLY CHANGES if (i==0 || i>=h+1 || j>=h) continue; peaks_prime[i] -=1; flats_prime[i-1]+=2; flats_prime[j] -=1; flats_prime[j+1]+=1; break; case 2: // FP(i,j) // TEST IF FEASIBLE AND APPLY CHANGES if (i==0 || j>=h) continue; peaks_prime[i] +=1; flats_prime[i-1]-=2; flats_prime[j] +=1; flats_prime[j+1]-=1; break; case 3: // FV(i,j) // TEST IF FEASIBLE AND APPLY CHANGES if (i==0 || i>=h+1 || j>=h) continue; flats_prime[i] -=2; peaks_prime[i] +=1; flats_prime[j] -=1; flats_prime[j+1]+=1; break; case 4: // VF(i,j) // TEST IF FEASIBLE AND APPLY CHANGES if (i==0 || i>=h+1 || j>=h) continue; flats_prime[i] +=2; peaks_prime[i] -=1; flats_prime[j] +=1; flats_prime[j+1]-=1; break; case 5: // FF(i,j) // TEST IF FEASIBLE AND APPLY CHANGES if (i==0 || i>=h+1 || j>=h) continue; flats_prime[i] -=1; flats_prime[i-1]+=1; flats_prime[j] -=1; flats_prime[j+1]+=1; break; case 6: // FF(i,j) // TEST IF FEASIBLE AND APPLY CHANGES if (i==0 || i>=h+1 || j>=h) continue; flats_prime[i] +=1; flats_prime[i-1]-=1; flats_prime[j] +=1; flats_prime[j+1]-=1; break; case 7: // PV(i,j) // TEST IF FEASIBLE AND APPLY CHANGES if (i==0 || i>=h+1 || j==0 || j>=h+1) continue; peaks_prime[i] -=1; flats_prime[i-1]+=2; peaks_prime[j] -=1; flats_prime[j] +=2; break; case 8: // VP(i,j) // TEST IF FEASIBLE AND APPLY CHANGES if (i==0 || j==0 || j>=h+1) continue; peaks_prime[i] +=1; flats_prime[i-1]-=2; peaks_prime[j] +=1; flats_prime[j] -=2; break; } // TEST IF THE STEP WAS FEASIBLE (NEW SEQUENCE IS VALID) bool feasible=true; for (int k=0; k<=h+1; k++) { if (flats_prime[k]<0 || peaks_prime[k]<0) feasible=false; } if (peaks_prime[0]!=0) feasible=false; for (int k=1; k0) feasible=false; if (peaks_prime[h+1]==0 && flats_prime[h+1]>0) feasible=false; if (peaks_prime[h]==0 && peaks_prime[h+1]>0) feasible=false; if (!feasible) continue; // COMPUTE THE TRANSITION PROBABILITY 'P = m(a')/m(a)' double P = 1.0; P *= binomial(peaks_prime[1]+flats_prime[0],flats_prime[0]); P /= binomial(peaks[1]+flats[0],flats[0]); for (int k=1; k<=h; k++) { P *= binomial(peaks_prime[k+1]+flats_prime[k]+peaks_prime[k]-1,peaks_prime[k]-1); P /= binomial(peaks[k+1]+flats[k]+peaks[k]-1,peaks[k]-1); P *= binomial(peaks_prime[k+1]+flats_prime[k],flats_prime[k]); P /= binomial(peaks[k+1]+flats[k],flats[k]); } P *= binomial(flats_prime[h+1]+peaks_prime[h+1]-1,peaks_prime[h+1]-1); P /= binomial(flats[h+1]+peaks[h+1]-1,peaks[h+1]-1); // COMPUTE THE TRANSITION PROBABILITY 'P *= perm(a')/perm(a)' P *= prob_perm(flats, flats_prime, peaks, peaks_prime); // ACCEPT OR REJECT THE MOVE if (randprob(gen)>min(1.0,P)) continue; if (peaks_prime[h+1]>0) { peaks_prime.push_back(0); flats_prime.push_back(0); h++; } else if (peaks_prime[h]==0) { peaks_prime.pop_back(); flats_prime.pop_back(); h--; } flats = flats_prime; peaks = peaks_prime; } // RETURN BUILDING BLOCK SEQUENCE for (int i=0; i<=h; i++) { if (i==0) { cout << flats[i];} else { cout << " " << peaks[i] << " " << flats[i];} } cout << endl; }