234 lines
6.4 KiB
C#
Raw Normal View History

// Copyright (C) 2016 Maxim Gumin, The MIT License (MIT)
2016-09-30 14:54:00 +04:00
using System;
abstract class Model
{
2019-09-27 12:58:58 +03:00
protected bool[][] wave;
protected int[][][] propagator;
int[][][] compatible;
protected int[] observed;
(int, int)[] stack;
int stacksize, observedSoFar;
2019-09-27 12:58:58 +03:00
protected int MX, MY, T, N;
protected bool periodic, ground;
2019-09-27 12:58:58 +03:00
protected double[] weights;
double[] weightLogWeights, distribution;
2019-09-27 12:58:58 +03:00
protected int[] sumsOfOnes;
2019-09-27 12:58:58 +03:00
double sumOfWeights, sumOfWeightLogWeights, startingEntropy;
protected double[] sumsOfWeights, sumsOfWeightLogWeights, entropies;
2019-09-27 12:58:58 +03:00
public enum Heuristic { Entropy, MRV, Scanline };
Heuristic heuristic;
protected Model(int width, int height, int N, bool periodic, Heuristic heuristic)
2019-09-27 12:58:58 +03:00
{
MX = width;
MY = height;
this.N = N;
this.periodic = periodic;
this.heuristic = heuristic;
2019-09-27 12:58:58 +03:00
}
void Init()
{
wave = new bool[MX * MY][];
2019-09-27 12:58:58 +03:00
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];
}
distribution = new double[T];
observed = new int[MX * MY];
2019-09-27 12:58:58 +03:00
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;
sumsOfOnes = new int[MX * MY];
sumsOfWeights = new double[MX * MY];
sumsOfWeightLogWeights = new double[MX * MY];
entropies = new double[MX * MY];
2019-09-27 12:58:58 +03:00
stack = new (int, int)[wave.Length * T];
stacksize = 0;
}
public bool Run(int seed, int limit)
2019-09-27 12:58:58 +03:00
{
if (wave == null) Init();
2019-09-27 12:58:58 +03:00
Clear();
Random random = new(seed);
for (int l = 0; l < limit || limit < 0; l++)
2019-09-27 12:58:58 +03:00
{
int node = NextUnobservedNode(random);
if (node >= 0)
{
Observe(node, random);
bool success = Propagate();
if (!success) return false;
}
else
{
for (int i = 0; i < wave.Length; i++) for (int t = 0; t < T; t++) if (wave[i][t]) { observed[i] = t; break; }
return true;
}
}
2019-09-27 12:58:58 +03:00
return true;
}
2019-09-27 12:58:58 +03:00
int NextUnobservedNode(Random random)
{
if (heuristic == Heuristic.Scanline)
{
for (int i = observedSoFar; i < wave.Length; i++)
{
if (!periodic && (i % MX + N > MX || i / MX + N > MY)) continue;
if (sumsOfOnes[i] > 1)
{
observedSoFar = i + 1;
return i;
}
}
return -1;
}
2021-10-21 13:23:03 +03:00
double min = 1E+4;
int argmin = -1;
for (int i = 0; i < wave.Length; i++)
{
if (!periodic && (i % MX + N > MX || i / MX + N > MY)) continue;
int remainingValues = sumsOfOnes[i];
double entropy = heuristic == Heuristic.Entropy ? entropies[i] : remainingValues;
if (remainingValues > 1 && entropy <= min)
2019-09-27 12:58:58 +03:00
{
double noise = 1E-6 * random.NextDouble();
if (entropy + noise < min)
{
min = entropy + noise;
argmin = i;
}
}
}
return argmin;
}
2019-09-27 12:58:58 +03:00
void Observe(int node, Random random)
{
bool[] w = wave[node];
for (int t = 0; t < T; t++) distribution[t] = w[t] ? weights[t] : 0.0;
2019-09-27 12:58:58 +03:00
int r = distribution.Random(random.NextDouble());
for (int t = 0; t < T; t++) if (w[t] != (t == r)) Ban(node, t);
2019-09-27 12:58:58 +03:00
}
bool Propagate()
2019-09-27 12:58:58 +03:00
{
while (stacksize > 0)
{
(int i1, int t1) = stack[stacksize - 1];
2019-09-27 12:58:58 +03:00
stacksize--;
int x1 = i1 % MX;
int y1 = i1 / MX;
2019-09-27 12:58:58 +03:00
for (int d = 0; d < 4; d++)
{
int x2 = x1 + dx[d];
int y2 = y1 + dy[d];
if (!periodic && (x2 < 0 || y2 < 0 || x2 + N > MX || y2 + N > MY)) continue;
2019-09-27 12:58:58 +03:00
if (x2 < 0) x2 += MX;
else if (x2 >= MX) x2 -= MX;
if (y2 < 0) y2 += MY;
else if (y2 >= MY) y2 -= MY;
2019-09-27 12:58:58 +03:00
int i2 = x2 + y2 * MX;
int[] p = propagator[d][t1];
2019-09-27 12:58:58 +03:00
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);
}
}
}
return sumsOfOnes[0] > 0;
2019-09-27 12:58:58 +03:00
}
void Ban(int i, int t)
2019-09-27 12:58:58 +03:00
{
wave[i][t] = false;
int[] comp = compatible[i][t];
for (int d = 0; d < 4; d++) comp[d] = 0;
stack[stacksize] = (i, t);
stacksize++;
sumsOfOnes[i] -= 1;
sumsOfWeights[i] -= weights[t];
sumsOfWeightLogWeights[i] -= weightLogWeights[t];
double sum = sumsOfWeights[i];
entropies[i] = Math.Log(sum) - sumsOfWeightLogWeights[i] / sum;
}
void Clear()
2019-09-27 12:58:58 +03:00
{
for (int i = 0; i < wave.Length; i++)
{
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;
observed[i] = -1;
2019-09-27 12:58:58 +03:00
}
observedSoFar = 0;
if (ground)
{
for (int x = 0; x < MX; x++)
{
for (int t = 0; t < T - 1; t++) Ban(x + (MY - 1) * MX, t);
for (int y = 0; y < MY - 1; y++) Ban(x + y * MX, T - 1);
}
Propagate();
}
2019-09-27 12:58:58 +03:00
}
public abstract void Save(string filename);
2019-09-27 12:58:58 +03:00
protected static int[] dx = { -1, 0, 1, 0 };
protected static int[] dy = { 0, 1, 0, -1 };
2019-09-27 12:58:58 +03:00
static int[] opposite = { 2, 3, 0, 1 };
2016-09-30 14:54:00 +04:00
}