From e2302deb9f435d0308d06780138eeff2659c81e4 Mon Sep 17 00:00:00 2001
From: Vincent Beffara <vbeffara@gmail.com>
Date: Wed, 11 May 2016 15:36:48 +0200
Subject: [PATCH] Move the domino shuffle into the 2D folder where it belongs

---
 tests/shuffle.cpp | 146 ----------------------------------------------
 1 file changed, 146 deletions(-)
 delete mode 100644 tests/shuffle.cpp

diff --git a/tests/shuffle.cpp b/tests/shuffle.cpp
deleted file mode 100644
index 3dd3b072..00000000
--- a/tests/shuffle.cpp
+++ /dev/null
@@ -1,146 +0,0 @@
-// Initial code in Python by Sunil Chhita, ported to C++ by VB.
-
-#include <vb/Array.h>
-#include <vb/Hub.h>
-#include <fstream>
-
-using namespace vb; using namespace std;
-
-auto weights (int m, const Array<double> & p) {
-	int n = m*lcm(p.ww,p.hh);
-	Array<double> A (n,n); for (auto z : coos(A)) A[z] = p.atp(z); return A;
-}
-
-auto dupe (const Array<double> & a) { Array<double> aa (2*a.ww,2*a.hh); for (auto z : coos(aa)) aa[z] = a[z/2]; return aa; }
-
-auto twoperiodic  	(int m, double a, double b)          	{ return weights (m, dupe (Array<double> ({{a,b},{b,a}}))); }
-auto threeperiodic	(int m, double a, double b, double c)	{ return weights (m, dupe (Array<double> ({{a,b,c},{b,c,a},{c,a,b}}))); }
-
-auto threebytwo(int m, double b) { Array<double> p (6,6,1); p[coo(0,0)]=b; p[coo(4,2)]=b; p[coo(2,4)]=b; return weights (m,p); }
-
-auto d3p (const Array<double> & x1) {
-	int n = x1.ww; vector<Array<pair<double,double>>> A (n/2);
-
-	A[0].resize(n,n); for (auto z : coos(x1)) if (double w = x1[z]) A[0][z] = {w,0}; else A[0][z] = {1,1};
-
-	for (int k=0; k<n/2-1; ++k) {
-		A[k+1].resize(n-2*k-2,n-2*k-2); for (auto z : coos(A[k+1])) {
-			int i=z.x, j=z.y, ii = i+2*(i%2), jj = j+2*(j%2);
-			double a20, a21;
-			if (A[k][coo(ii,jj)].second + A[k][coo(i+1,j+1)].second == A[k][coo(ii,j+1)].second + A[k][coo(i+1,jj)].second) {
-				a20 = A[k][coo(ii,jj)].first * A[k][coo(i+1,j+1)].first + A[k][coo(ii,j+1)].first * A[k][coo(i+1,jj)].first;
-				a21 = A[k][coo(ii,jj)].second + A[k][coo(i+1,j+1)].second;
-			} else if (A[k][coo(ii,jj)].second + A[k][coo(i+1,j+1)].second < A[k][coo(ii,j+1)].second + A[k][coo(i+1,jj)].second) {
-				a20 = A[k][coo(ii,jj)].first * A[k][coo(i+1,j+1)].first;
-				a21 = A[k][coo(ii,jj)].second + A[k][coo(i+1,j+1)].second;
-			} else {
-				a20 = A[k][coo(ii,j+1)].first * A[k][coo(i+1,jj)].first;
-				a21 = A[k][coo(ii,j+1)].second + A[k][coo(i+1,jj)].second;
-			}
-			auto & a1 = A[k][coo(ii,jj)]; A[k+1][z] = { a1.first/a20, a1.second-a21 };
-		}
-	}
-	return A;
-}
-
-auto probs (const Array<double> & x1) {
-	auto a0 = d3p (x1); int n = a0.size(); vector<Array<double>> A (n);
-
-	for (int k=0; k<n; ++k) {
-		A[k].resize (k+1,k+1); for (auto z : coos(A[k])) {
-			int i=z.x, j=z.y;
-			if (a0[n-k-1][coo(2*i,2*j)].second + a0[n-k-1][coo(2*i+1,2*j+1)].second > a0[n-k-1][coo(2*i+1,2*j)].second + a0[n-k-1][coo(2*i,2*j+1)].second)
-				A[k][z] = 0;
-			else if (a0[n-k-1][coo(2*i,2*j)].second + a0[n-k-1][coo(2*i+1,2*j+1)].second < a0[n-k-1][coo(2*i+1,2*j)].second + a0[n-k-1][coo(2*i,2*j+1)].second)
-				A[k][z] = 1;
-			else
-				A[k][z] = (a0[n-k-1][coo(2*i+1,2*j+1)].first * a0[n-k-1][coo(2*i,2*j)].first /
-					(a0[n-k-1][coo(2*i+1,2*j+1)].first * a0[n-k-1][coo(2*i,2*j)].first + a0[n-k-1][coo(2*i+1,2*j)].first*a0[n-k-1][coo(2*i,2*j+1)].first));
-		}
-	}
-	return A;
-}
-
-auto delslide (const Array<int> & x1) {
-	int n = x1.ww; Array<int> a0 (n+2,n+2);
-	for (int i=0; i<n+2; ++i)
-		for (int j=0; j<n+2; ++j)
-			a0[coo(i,j)] = ((i==0) || (i==n+1) || (j==0) || (j==n+1)) ? 0 : x1[coo(i-1,j-1)];
-	for (int i=0; i<n/2; ++i) {
-		for (int j=0; j<n/2; ++j) {
-			if ((a0[coo(2*i,2*j)]==1) && (a0[coo(2*i+1,2*j+1)]==1)) { a0[coo(2*i,2*j)]=0; a0[coo(2*i+1,2*j+1)]=0; }
-			else if ((a0[coo(2*i,2*j+1)]==1) && (a0[coo(2*i+1,2*j)]==1)) { a0[coo(2*i+1,2*j)]=0; a0[coo(2*i,2*j+1)]=0; }
-		}
-	}
-	for (int i=0; i<n/2+1; ++i) {
-		for (int j=0; j<n/2+1; ++j) {
-			if (a0[coo(2*i+1,2*j+1)]==1) { a0[coo(2*i,2*j)]=1; a0[coo(2*i+1,2*j+1)]=0; }
-			else if (a0[coo(2*i,2*j)]==1) { a0[coo(2*i,2*j)]=0; a0[coo(2*i+1,2*j+1)]=1; }
-			else if (a0[coo(2*i+1,2*j)]==1) { a0[coo(2*i,2*j+1)]=1; a0[coo(2*i+1,2*j)]=0; }
-			else if (a0[coo(2*i,2*j+1)]==1) { a0[coo(2*i+1,2*j)]=1; a0[coo(2*i,2*j+1)]=0; }
-		}
-	}
-	return a0;
-}
-
-auto create (Array<int> & x0, const Array<double> & p) {
-	int n = x0.ww;
-	for (int i=0; i<n/2; ++i) {
-		for (int j=0; j<n/2; ++j) {
-			if ((x0[coo(2*i,2*j)]==0) && (x0[coo(2*i+1,2*j)]==0) && (x0[coo(2*i,2*j+1)]==0) && (x0[coo(2*i+1,2*j+1)]==0)) {
-				bool a1,a2,a3,a4;
-				if (j>0) a1 = (x0[coo(2*i,2*j-1)]==0) && (x0[coo(2*i+1,2*j-1)]==0); else a1 = true;
-				if (j<n/2-1) a2 = (x0[coo(2*i,2*j+2)]==0) && (x0[coo(2*i+1,2*j+2)]==0); else a2 = true;
-				if (i>0) a3 = (x0[coo(2*i-1,2*j)]==0) && (x0[coo(2*i-1,2*j+1)]==0); else a3 = true;
-				if (i<n/2-1) a4 = (x0[coo(2*i+2,2*j)]==0) && (x0[coo(2*i+2,2*j+1)]==0); else a4 = true;
-				if (a1 && a2 && a3 && a4) {
-					if (prng.bernoulli(p[coo(i,j)])) { x0[coo(2*i,2*j)]=1; x0[coo(2*i+1,2*j+1)]=1; }
-					else { x0[coo(2*i+1,2*j)]=1; x0[coo(2*i,2*j+1)]=1; }
-				}
-			}
-		}
-	}
-}
-
-auto aztecgen (const Array<double> & xr) {
-	auto x0 = probs(xr); int n = x0.size(); Array<int> a1;
-	if (prng.bernoulli(x0[0][coo(0,0)])) a1 = Array<int> ({{1,0},{0,1}}); else a1 = Array<int> ({{0,1},{1,0}});
-	for (int i=0; i<n-1; ++i) { a1=delslide(a1); create(a1,x0[i+1]); }
-	return a1;
-}
-
-auto height (const Array<int> & A) {
-	int m = A.ww/2; Array<int> h (m+1,m+1,0); int z=0;
-	for (int x=0; x<m; ++x) { z += 4*A[coo(2*x,0)] + 4*A[coo(2*x+1,0)] - 2; h[coo(x+1,0)] = z; }
-	for (int y=0; y<m; ++y) {
-		int z = h[coo(0,y)] + 4*A[coo(0,2*y)] + 4*A[coo(0,2*y+1)] - 2; h[coo(0,y+1)] = z;
-		for (int x=0; x<m; ++x) { z -= 4*A[coo(2*x,2*y+1)] + 4*A[coo(2*x+1,2*y+1)] - 2; h[coo(x+1,y+1)] = z; }
-	}
-	return h;
-}
-
-int main (int argc, char ** argv) {
-	H.init ("Domino shuffle", argc,argv, "m=50,a=1,b=.5,c=.3,w=unif");
-	int m = H['m'];
-
-	Array<double> TP;
-	if (H['w'] == "unif") TP = weights (m, Array<double> (1,1,1)); else
-	if (H['w'] == "two") TP = twoperiodic (m,H['a'],H['b']); else
-	if (H['w'] == "three") TP = threeperiodic (m,H['a'],H['b'],H['c']); else
-	if (H['w'] == "kenyon") TP = threebytwo (m,H['b']); else {
-		cerr << "No such weight, \"" << string(H['w']) << "\".\n"; exit(1);
-	}
-
-	auto A1 = aztecgen(TP); auto H1 = height(A1);
-
-	string name = H.dir + H.title; ofstream asy (name + ".asy");
-	for (auto z : coos(A1)) if (A1[z]) {
-		double eps = ((z.x+z.y)%2) ? .5 : -.5;
-		asy << "draw ((" << z.x-.5 << "," << z.y-eps << ")--(" << z.x+.5 << "," << z.y+eps << "), gray (" << TP[z]/1.3 << "));" << endl;
-	}
-
-	ofstream dat (name + ".dat");
-	for (auto z : coos(H1)) { dat << H1[z] << " "; if (z.x == H1.ww-1) dat << endl; }
-
-	system (("asy -fpdf -o " + H.dir + " \"" + name + ".asy\"").c_str());
-}