From fad1066b5000f7e9fbda0ef81bbea56799686670 Mon Sep 17 00:00:00 2001 From: ExUtumno Date: Tue, 27 Mar 2018 00:29:55 +0300 Subject: [PATCH] optimizations of Fehr and Courant, x20 speed-up for the overlapping model --- Main.cs | 5 ++ Model.cs | 175 +++++++++++++++++++++++++++++++------------- OverlappingModel.cs | 86 ++++------------------ SimpleTiledModel.cs | 89 ++-------------------- 4 files changed, 150 insertions(+), 205 deletions(-) diff --git a/Main.cs b/Main.cs index abd5415..cd1748f 100644 --- a/Main.cs +++ b/Main.cs @@ -8,11 +8,14 @@ The above copyright notice and this permission notice shall be included in all c using System; using System.Xml.Linq; +using System.Diagnostics; static class Program { static void Main() { + Stopwatch sw = Stopwatch.StartNew(); + Random random = new Random(); XDocument xdoc = XDocument.Load("samples.xml"); @@ -52,5 +55,7 @@ static void Main() counter++; } + + Console.WriteLine($"time = {sw.ElapsedMilliseconds}"); } } diff --git a/Model.cs b/Model.cs index 2c12fdb..4d2c4dd 100644 --- a/Model.cs +++ b/Model.cs @@ -11,34 +11,64 @@ The above copyright notice and this permission notice shall be included in all c abstract class Model { protected bool[][] wave; - protected double[] stationary; + + protected int[][][] propagator; + int[][][] compatible; protected int[] observed; - protected bool[] changes; - protected int[] stack; - protected int stacksize; + Tuple[] stack; + int stacksize; protected Random random; protected int FMX, FMY, T; protected bool periodic; - double[] logProb; - double logT; + protected double[] weights; + double[] weightLogWeights; + + int[] sumsOfOnes; + double sumOfWeights, sumOfWeightLogWeights, startingEntropy; + double[] sumsOfWeights, sumsOfWeightLogWeights, entropies; - protected Model(int width, int height) + protected Model(int width, int height) { FMX = width; FMY = height; + } + void Init() + { wave = new bool[FMX * FMY][]; - changes = new bool[FMX * FMY]; + compatible = new int[wave.Length][][]; + for (int i = 0; i < wave.Length; i++) + { + wave[i] = new bool[T]; + compatible[i] = new int[T][]; + for (int t = 0; t < T; t++) compatible[i][t] = new int[4]; + } + + weightLogWeights = new double[T]; + sumOfWeights = 0; + sumOfWeightLogWeights = 0; + + for (int t = 0; t < T; t++) + { + weightLogWeights[t] = weights[t] * Math.Log(weights[t]); + sumOfWeights += weights[t]; + sumOfWeightLogWeights += weightLogWeights[t]; + } + + startingEntropy = Math.Log(sumOfWeights) - sumOfWeightLogWeights / sumOfWeights; - stack = new int[FMX * FMY]; + sumsOfOnes = new int[FMX * FMY]; + sumsOfWeights = new double[FMX * FMY]; + sumsOfWeightLogWeights = new double[FMX * FMY]; + entropies = new double[FMX * FMY]; + + stack = new Tuple[wave.Length * T]; stacksize = 0; } - protected abstract void Propagate(); - bool? Observe() { double min = 1E+3; @@ -46,37 +76,20 @@ protected Model(int width, int height) for (int i = 0; i < wave.Length; i++) { - if (OnBoundary(i)) continue; + if (OnBoundary(i % FMX, i / FMX)) continue; - bool[] w = wave[i]; - int amount = 0; - double sum = 0; + int amount = sumsOfOnes[i]; + if (amount == 0) return false; - for (int t = 0; t < T; t++) if (w[t]) + double entropy = entropies[i]; + if (amount > 1 && entropy <= min) + { + double noise = 1E-6 * random.NextDouble(); + if (entropy + noise < min) { - amount += 1; - sum += stationary[t]; + min = entropy + noise; + argmin = i; } - - if (sum == 0) return false; - - double noise = 1E-6 * random.NextDouble(); - - double entropy; - if (amount == 1) entropy = 0; - else if (amount == T) entropy = logT; - else - { - double mainSum = 0; - double logSum = Math.Log(sum); - for (int t = 0; t < T; t++) if (w[t]) mainSum += stationary[t] * logProb[t]; - entropy = logSum - mainSum / sum; - } - - if (entropy > 0 && entropy + noise < min) - { - min = entropy + noise; - argmin = i; } } @@ -88,19 +101,56 @@ protected Model(int width, int height) } double[] distribution = new double[T]; - for (int t = 0; t < T; t++) distribution[t] = wave[argmin][t] ? stationary[t] : 0; + for (int t = 0; t < T; t++) distribution[t] = wave[argmin][t] ? weights[t] : 0; int r = distribution.Random(random.NextDouble()); - for (int t = 0; t < T; t++) wave[argmin][t] = t == r; - Change(argmin); + + bool[] w = wave[argmin]; + for (int t = 0; t < T; t++) if (w[t] != (t == r)) Ban(argmin, t); return null; } + protected void Propagate() + { + while (stacksize > 0) + { + var e1 = stack[stacksize - 1]; + stacksize--; + + int i1 = e1.Item1; + int x1 = i1 % FMX, y1 = i1 / FMX; + bool[] w1 = wave[i1]; + + for (int d = 0; d < 4; d++) + { + int dx = DX[d], dy = DY[d]; + int x2 = x1 + dx, y2 = y1 + dy; + if (OnBoundary(x2, y2)) continue; + + if (x2 < 0) x2 += FMX; + else if (x2 >= FMX) x2 -= FMX; + if (y2 < 0) y2 += FMY; + else if (y2 >= FMY) y2 -= FMY; + + int i2 = x2 + y2 * FMX; + int[] p = propagator[d][e1.Item2]; + int[][] compat = compatible[i2]; + + for (int l = 0; l < p.Length; l++) + { + int t2 = p[l]; + int[] comp = compat[t2]; + + comp[d]--; + if (comp[d] == 0) Ban(i2, t2); + } + } + } + } + public bool Run(int seed, int limit) { - logT = Math.Log(T); - logProb = new double[T]; - for (int t = 0; t < T; t++) logProb[t] = Math.Log(stationary[t]); + if (wave == null) Init(); Clear(); random = new Random(seed); @@ -115,24 +165,47 @@ public bool Run(int seed, int limit) return true; } - protected void Change(int i) + protected void Ban(int i, int t) { - if (changes[i]) return; + wave[i][t] = false; - stack[stacksize] = i; + int[] comp = compatible[i][t]; + for (int d = 0; d < 4; d++) comp[d] = 0; + stack[stacksize] = new Tuple(i, t); stacksize++; - changes[i] = true; + + double sum = sumsOfWeights[i]; + entropies[i] += sumsOfWeightLogWeights[i] / sum - Math.Log(sum); + + sumsOfOnes[i] -= 1; + sumsOfWeights[i] -= weights[t]; + sumsOfWeightLogWeights[i] -= weightLogWeights[t]; + + sum = sumsOfWeights[i]; + entropies[i] -= sumsOfWeightLogWeights[i] / sum - Math.Log(sum); } protected virtual void Clear() { for (int i = 0; i < wave.Length; i++) { - for (int t = 0; t < T; t++) wave[i][t] = true; - changes[i] = false; + for (int t = 0; t < T; t++) + { + wave[i][t] = true; + for (int d = 0; d < 4; d++) compatible[i][t][d] = propagator[opposite[d]][t].Length; + } + + sumsOfOnes[i] = weights.Length; + sumsOfWeights[i] = sumOfWeights; + sumsOfWeightLogWeights[i] = sumOfWeightLogWeights; + entropies[i] = startingEntropy; } } - protected abstract bool OnBoundary(int i); + protected abstract bool OnBoundary(int x, int y); public abstract System.Drawing.Bitmap Graphics(); + + protected static int[] DX = { -1, 0, 1, 0 }; + protected static int[] DY = { 0, 1, 0, -1 }; + static int[] opposite = { 2, 3, 0, 1 }; } diff --git a/OverlappingModel.cs b/OverlappingModel.cs index 70c82d2..ce13fd0 100644 --- a/OverlappingModel.cs +++ b/OverlappingModel.cs @@ -13,9 +13,7 @@ The above copyright notice and this permission notice shall be included in all c class OverlappingModel : Model { - int[][][][] propagator; int N; - byte[][] patterns; List colors; int ground; @@ -123,21 +121,17 @@ byte[] patternFromIndex(long ind) T = weights.Count; this.ground = (ground + T) % T; - patterns = new byte[T][]; - stationary = new double[T]; - propagator = new int[2 * N - 1][][][]; + base.weights = new double[T]; int counter = 0; foreach (long w in ordering) { patterns[counter] = patternFromIndex(w); - stationary[counter] = weights[w]; + base.weights[counter] = weights[w]; counter++; } - for (int i = 0; i < wave.Length; i++) wave[i] = new bool[T]; - bool agrees(byte[] p1, byte[] p2, int dx, int dy) { int xmin = dx < 0 ? 0 : dx, xmax = dx < 0 ? dx + N : N, ymin = dy < 0 ? 0 : dy, ymax = dy < 0 ? dy + N : N; @@ -145,67 +139,21 @@ bool agrees(byte[] p1, byte[] p2, int dx, int dy) return true; }; - for (int x = 0; x < 2 * N - 1; x++) + propagator = new int[4][][]; + for (int d = 0; d < 4; d++) { - propagator[x] = new int[2 * N - 1][][]; - for (int y = 0; y < 2 * N - 1; y++) + propagator[d] = new int[T][]; + for (int t = 0; t < T; t++) { - propagator[x][y] = new int[T][]; - for (int t = 0; t < T; t++) - { - List list = new List(); - for (int t2 = 0; t2 < T; t2++) if (agrees(patterns[t], patterns[t2], x - N + 1, y - N + 1)) list.Add(t2); - propagator[x][y][t] = new int[list.Count]; - for (int c = 0; c < list.Count; c++) propagator[x][y][t][c] = list[c]; - } + List list = new List(); + for (int t2 = 0; t2 < T; t2++) if (agrees(patterns[t], patterns[t2], DX[d], DY[d])) list.Add(t2); + propagator[d][t] = new int[list.Count]; + for (int c = 0; c < list.Count; c++) propagator[d][t][c] = list[c]; } } } - protected override bool OnBoundary(int i) => !periodic && (i % FMX + N > FMX || i / FMX + N > FMY); - - override protected void Propagate() - { - while (stacksize > 0) - { - int i1 = stack[stacksize - 1]; - stacksize--; - changes[i1] = false; - - bool[] w1 = wave[i1]; - int x1 = i1 % FMX, y1 = i1 / FMX; - - for (int dx = -N + 1; dx < N; dx++) for (int dy = -N + 1; dy < N; dy++) - { - int x2 = x1 + dx; - if (x2 < 0) x2 += FMX; - else if (x2 >= FMX) x2 -= FMX; - - int y2 = y1 + dy; - if (y2 < 0) y2 += FMY; - else if (y2 >= FMY) y2 -= FMY; - - if (!periodic && (x2 + N > FMX || y2 + N > FMY)) continue; - - int i2 = x2 + y2 * FMX; - bool[] w2 = wave[i2]; - int[][] prop = propagator[N - 1 - dx][N - 1 - dy]; - - for (int t2 = 0; t2 < T; t2++) if (w2[t2]) - { - bool b = false; - int[] p = prop[t2]; - for (int l = 0; l < p.Length && !b; l++) b = w1[p[l]]; - - if (!b) - { - Change(i2); - w2[t2] = false; - } - } - } - } - } + protected override bool OnBoundary(int x, int y) => !periodic && (x + N > FMX || y + N > FMY || x < 0 || y < 0); public override Bitmap Graphics() { @@ -241,7 +189,7 @@ public override Bitmap Graphics() if (sy < 0) sy += FMY; int s = sx + sy * FMX; - if (OnBoundary(s)) continue; + if (OnBoundary(sx, sy)) continue; for (int t = 0; t < T; t++) if (wave[s][t]) { contributors++; @@ -271,14 +219,8 @@ protected override void Clear() { for (int x = 0; x < FMX; x++) { - for (int t = 0; t < T; t++) if (t != ground) wave[x + (FMY - 1) * FMX][t] = false; - Change(x + (FMY - 1) * FMX); - - for (int y = 0; y < FMY - 1; y++) - { - wave[x + y * FMX][ground] = false; - Change(x + y * FMX); - } + for (int t = 0; t < T; t++) if (t != ground) Ban(x + (FMY - 1) * FMX, t); + for (int y = 0; y < FMY - 1; y++) Ban(x + y * FMX, ground); } Propagate(); diff --git a/SimpleTiledModel.cs b/SimpleTiledModel.cs index a07a471..1ea74fc 100644 --- a/SimpleTiledModel.cs +++ b/SimpleTiledModel.cs @@ -15,8 +15,6 @@ The above copyright notice and this permission notice shall be included in all c class SimpleTiledModel : Model { - int[][][] propagator; - List tiles; List tilenames; int tilesize; @@ -143,10 +141,10 @@ Color[] tile (Func f) } T = action.Count; - stationary = tempStationary.ToArray(); + weights = tempStationary.ToArray(); - var tempPropagator = new bool[4][][]; propagator = new int[4][][]; + var tempPropagator = new bool[4][][]; for (int d = 0; d < 4; d++) { tempPropagator[d] = new bool[T][]; @@ -154,8 +152,6 @@ Color[] tile (Func f) for (int t = 0; t < T; t++) tempPropagator[d][t] = new bool[T]; } - for (int i = 0; i < wave.Length; i++) wave[i] = new bool[T]; - foreach (XElement xneighbor in xroot.Element("neighbors").Elements("neighbor")) { string[] left = xneighbor.Get("left").Split(new char[] { ' ' }, StringSplitOptions.RemoveEmptyEntries); @@ -203,78 +199,7 @@ Color[] tile (Func f) } } - protected override void Propagate() - { - while (stacksize > 0) - { - int i1 = stack[stacksize - 1]; - changes[i1] = false; - stacksize--; - - bool[] w1 = wave[i1]; - int x1 = i1 % FMX, y1 = i1 / FMX; - - for (int d = 0; d < 4; d++) - { - int x2 = x1, y2 = y1; - if (d == 0) - { - if (x1 == FMX - 1) - { - if (!periodic) continue; - else x2 = 0; - } - else x2 = x1 + 1; - } - else if (d == 1) - { - if (y1 == 0) - { - if (!periodic) continue; - else y2 = FMY - 1; - } - else y2 = y1 - 1; - } - else if (d == 2) - { - if (x1 == 0) - { - if (!periodic) continue; - else x2 = FMX - 1; - } - else x2 = x1 - 1; - } - else - { - if (y1 == FMY - 1) - { - if (!periodic) continue; - else y2 = 0; - } - else y2 = y1 + 1; - } - - int i2 = x2 + y2 * FMX; - bool[] w2 = wave[i2]; - int[][] prop = propagator[d]; - - for (int t2 = 0; t2 < T; t2++) if (w2[t2]) - { - bool b = false; - int[] p = prop[t2]; - for (int l = 0; l < p.Length && !b; l++) b = w1[p[l]]; - - if (!b) - { - Change(i2); - w2[t2] = false; - } - } - } - } - } - - protected override bool OnBoundary(int i) => false; + protected override bool OnBoundary(int x, int y) => !periodic && (x < 0 || y < 0 || x >= FMX || y >= FMY); public override Bitmap Graphics() { @@ -300,7 +225,7 @@ public override Bitmap Graphics() { bool[] a = wave[x + y * FMX]; int amount = (from b in a where b select 1).Sum(); - double lambda = 1.0 / (from t in Enumerable.Range(0, T) where a[t] select stationary[t]).Sum(); + double lambda = 1.0 / (from t in Enumerable.Range(0, T) where a[t] select weights[t]).Sum(); for (int yt = 0; yt < tilesize; yt++) for (int xt = 0; xt < tilesize; xt++) { @@ -311,9 +236,9 @@ public override Bitmap Graphics() for (int t = 0; t < T; t++) if (wave[x + y * FMX][t]) { Color c = tiles[t][xt + yt * tilesize]; - r += (double)c.R * stationary[t] * lambda; - g += (double)c.G * stationary[t] * lambda; - b += (double)c.B * stationary[t] * lambda; + r += (double)c.R * weights[t] * lambda; + g += (double)c.G * weights[t] * lambda; + b += (double)c.B * weights[t] * lambda; } bitmapData[x * tilesize + xt + (y * tilesize + yt) * FMX * tilesize] =