using System;
using System.Collections.Generic;
using System.IO;
using System.Linq;

namespace Amphora.RandomComparisons
{
    /// <summary>
    /// Checks the checks.  The statistics are compared with values computed independently by scipy, fractions and a
    /// separate Python implementation.  The generators and controls are run through a quick battery.  Every control
    /// must fail what it is known to fail, and no real generator may fail anything.
    /// </summary>
    public static class SelfCheck
    {
        private sealed class Reporter
        {
            private readonly TextWriter output;

            public int Checks;
            public int Failures;

            public Reporter(TextWriter output)
            {
                this.output = output;
            }

            public void Check(string name, bool passed, string detail = "")
            {
                Checks++;
                if (!passed)
                {
                    Failures++;
                    output.WriteLine("FAILED  " + name + (detail == "" ? "" : ": " + detail));
                }
            }

            public void Near(string name, double expected, double actual, double relativeTolerance)
            {
                double allowed = Math.Max(Math.Abs(expected) * relativeTolerance, 1e-300);
                Check(name, Math.Abs(expected - actual) <= allowed, "expected " + expected + " but was " + actual);
            }
        }

        private static readonly (double Statistic, double Degrees, double P)[] ChiSquareCases =
        {
            (65871.5, 65535, 0.3525158799117868),
            (16000.0, 16383, 0.03326359930771983),
            (262228.0, 262143, 0.9058307485331057),
            (50.0, 100, 1.3906610495232403e-05),
            (5.3, 6, 0.9880980009419972),
            (1060.3, 999, 0.17409941759607617),
            (4981699.5, 262143, 0.0),
            (3.0, 3, 0.7832503525421755),
            (0.5, 3, 0.16221717669064836),
            (12.0, 3, 0.014766321010719539),
            (66150.9, 65535, 0.08955455299533015),
            (64611.4, 65535, 0.010426067734684415),
            (16384.0, 16383, 0.9926537570520125),
            (1000000.0, 999999, 0.9990596839762254),
        };

        private static readonly (double Z, double P)[] NormalCases =
        {
            (0.0, 1.0),
            (0.5, 0.6170750774519738),
            (2.0, 0.0455002638963584),
            (3.0, 0.002699796063260186),
            (5.0, 5.733031437583867e-07),
            (8.0, 1.244192114854348e-15),
        };

        private static byte[] BitsOf(ulong seed, int length, Window window)
        {
            var generator = new SplitMix64(seed);
            var returnMe = new byte[length];
            for (int i = 0; i < length; i++)
            {
                returnMe[i] = (byte)Windows.EdgeBit(generator.NextUInt64(), window);
            }

            return returnMe;
        }

        private static ulong[] RowsOf(ulong seed, int size, Window window)
        {
            var generator = new SplitMix64(seed);
            var returnMe = new ulong[size];
            for (int i = 0; i < size; i++)
            {
                ulong word = generator.NextUInt64();
                returnMe[i] = size == 64 ? word : Windows.Take(word, window, size);
            }

            return returnMe;
        }

        private static double Mean(double[] probabilities)
        {
            double returnMe = 0;
            for (int total = 0; total < probabilities.Length; total++)
            {
                returnMe += total * probabilities[total];
            }

            return returnMe;
        }

        private static void CheckDiceDistributions(Reporter reporter)
        {
            double[] threeD6 = DistributionChecks.SumProbabilities(3, 6);
            reporter.Near("3d6 at 10", 0.125, threeD6[10], 1e-12);
            reporter.Near("3d6 at 3", 1.0 / 216, threeD6[3], 1e-12);
            reporter.Near("3d6 adds to one", 1, threeD6.Sum(), 1e-12);
            double[] dropLowest = DistributionChecks.KeptProbabilities(4, 6, 3, highest: true);
            reporter.Near("4d6 dropping the lowest at 18", 7.0 / 432, dropLowest[18], 1e-12);
            reporter.Near("4d6 dropping the lowest at 3", 1.0 / 1296, dropLowest[3], 1e-12);
            reporter.Near("4d6 dropping the lowest, mean", 12.244598765432098, Mean(dropLowest), 1e-12);
            reporter.Near("4d6 dropping the lowest adds to one", 1, dropLowest.Sum(), 1e-12);
            double[] advantage = DistributionChecks.KeptProbabilities(2, 20, 1, highest: true);
            reporter.Near("best of two d20 at 20", 0.0975, advantage[20], 1e-12);
            reporter.Near("best of two d20 at 1", 0.0025, advantage[1], 1e-12);
            double[] disadvantage = DistributionChecks.KeptProbabilities(3, 20, 1, highest: false);
            reporter.Near("worst of three d20 at 1", 0.142625, disadvantage[1], 1e-12);
            reporter.Near("worst of three d20 at 20", 0.000125, disadvantage[20], 1e-12);
            reporter.Near("worst of three d20, mean", 5.5125, Mean(disadvantage), 1e-12);
            double[] tenD6 = DistributionChecks.KeptProbabilities(10, 6, 5, highest: true);
            reporter.Near("best five of ten d6 at 30", 0.015461966703500482, tenD6[30], 1e-9);
            reporter.Near("best five of ten d6, mean", 24.179310264303798, Mean(tenD6), 1e-9);
            reporter.Near("best five of ten d6 adds to one", 1, tenD6.Sum(), 1e-9);
        }

        private static void CheckStatistics(Reporter reporter)
        {
            foreach (var testCase in ChiSquareCases)
            {
                reporter.Near("chi-square " + testCase.Statistic + " on " + testCase.Degrees + " df", testCase.P,
                    Stats.ChiSquareTwoSided(testCase.Statistic, testCase.Degrees), 1e-6);
            }

            foreach (var testCase in NormalCases)
            {
                reporter.Near("normal " + testCase.Z, testCase.P, Stats.NormalTwoSided(testCase.Z), 1e-6);
            }

            double[] eight = Stats.PoissonProbabilities(8, 16);
            reporter.Near("Poisson(8) at 0", 0.00033546262790251185, eight[0], 1e-9);
            reporter.Near("Poisson(8) at 5", 0.09160366159257921, eight[5], 1e-9);
            reporter.Near("Poisson(8) at 8", 0.13958653195059664, eight[8], 1e-9);
            reporter.Near("Poisson(8) at 15", 0.009025979411121507, eight[15], 1e-9);
            double[] four = Stats.PoissonProbabilities(4, 16);
            reporter.Near("Poisson(4) at 1", 0.07326255555493673, four[1], 1e-9);
            reporter.Near("Poisson(4) at 8", 0.029770181304863145, four[8], 1e-9);
            reporter.Near("32x32 full rank", 0.28878809515384113, Stats.BinaryRankProbability(32, 32), 1e-9);
            reporter.Near("32x32 rank 31", 0.5775761901732048, Stats.BinaryRankProbability(32, 31), 1e-9);
            reporter.Near("32x32 rank 30", 0.12835026442316666, Stats.BinaryRankProbability(32, 30), 1e-9);
            reporter.Near("64x64 full rank", 0.2887880950866024, Stats.BinaryRankProbability(64, 64), 1e-9);
            reporter.Near("64x64 rank 63", 0.5775761901732048, Stats.BinaryRankProbability(64, 63), 1e-9);
            reporter.Near("64x64 rank 62", 0.1283502644829344, Stats.BinaryRankProbability(64, 62), 1e-9);
            reporter.Near("worst of 8 at 0.01", 1 - Math.Pow(0.99, 8), Stats.WorstOf(0.01, 8), 1e-12);
            reporter.Near("worst of 8 at 1e-12", 8e-12, Stats.WorstOf(1e-12, 8), 1e-9);
            reporter.Near("worst of 8 at 1", 1, Stats.WorstOf(1, 8), 1e-12);
            CheckDiceDistributions(reporter);
            double statistic;
            double degrees = Stats.ChiSquareOfCategories(new double[] { 10, 1, 1, 8 }, new double[] { 9, 2, 2, 7 }, 5,
                out statistic);
            reporter.Check("categories merge until each is large enough", degrees == 1,
                "degrees " + degrees + " but expected 1");
            reporter.Near("pooled statistic", 1.0 / 9 + 1.0 / 11, statistic, 1e-12);
        }

        private static void CheckLinearAlgebra(Reporter reporter)
        {
            var machine = new BerlekampMassey(2000);
            reporter.Check("complexity of 50 bits", machine.Complexity(BitsOf(1, 50, Window.High), 50) == 26);
            reporter.Check("complexity of 300 bits", machine.Complexity(BitsOf(2, 300, Window.High), 300) == 150);
            reporter.Check("complexity of 2000 bits", machine.Complexity(BitsOf(3, 2000, Window.High), 2000) == 1001);
            reporter.Check("complexity of another 2000 bits",
                machine.Complexity(BitsOf(4, 2000, Window.High), 2000) == 1000);
            reporter.Check("complexity of 2000 low bits",
                machine.Complexity(BitsOf(5, 2000, Window.Low), 2000) == 1000);

            var lfsr = new byte[400];
            int[] start = { 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0 };
            for (int i = 0; i < lfsr.Length; i++)
            {
                lfsr[i] = i < start.Length ? (byte)start[i] : (byte)(lfsr[i - 3] ^ lfsr[i - 17]);
            }

            reporter.Check("a degree 17 register has complexity 17", machine.Complexity(lfsr, lfsr.Length) == 17);
            reporter.Check("zeros have complexity 0", machine.Complexity(new byte[100], 100) == 0);
            var lateOne = new byte[100];
            lateOne[99] = 1;
            reporter.Check("a lone late one has full complexity", machine.Complexity(lateOne, 100) == 100);
            var alternating = new byte[100];
            for (int i = 0; i < alternating.Length; i++)
            {
                alternating[i] = (byte)(i % 2);
            }

            reporter.Check("alternating bits have complexity 2", machine.Complexity(alternating, 100) == 2);
            reporter.Check("a machine can be reused", machine.Complexity(lfsr, lfsr.Length) == 17);

            int[] lowRanks = { 31, 32, 31, 31, 31, 30 };
            for (int seed = 1; seed <= lowRanks.Length; seed++)
            {
                reporter.Check("rank of 32x32 matrix " + seed,
                    BinaryRank.Rank(RowsOf((ulong)seed, 32, Window.Low), 32) == lowRanks[seed - 1]);
            }

            for (int seed = 1; seed <= 3; seed++)
            {
                reporter.Check("rank of 64x64 matrix " + seed,
                    BinaryRank.Rank(RowsOf((ulong)seed, 64, Window.High), 64) == 63);
            }

            var identity = new ulong[32];
            for (int i = 0; i < identity.Length; i++)
            {
                identity[i] = 1UL << i;
            }

            reporter.Check("the identity has full rank", BinaryRank.Rank(identity, 32) == 32);
            reporter.Check("zeros have rank 0", BinaryRank.Rank(new ulong[32], 32) == 0);
            identity[5] = identity[4];
            reporter.Check("a repeated row costs a rank", BinaryRank.Rank(identity, 32) == 31);

            int[] expectedCategories = { 0, 1, 3, 4, 5, 6 };
            int[] complexities = { 997, 998, 1000, 1001, 1002, 1003 };
            for (int i = 0; i < complexities.Length; i++)
            {
                reporter.Check("complexity " + complexities[i] + " of 2000 falls in category " +
                    expectedCategories[i], LinearComplexity.Category(complexities[i], 2000) == expectedCategories[i]);
            }

            reporter.Near("complexity categories cover everything", 1, LinearComplexity.Probabilities.Sum(), 1e-12);
            reporter.Check("popcount of zero", Bits.PopCount(0) == 0);
            reporter.Check("popcount of all ones", Bits.PopCount(ulong.MaxValue) == 64);
            reporter.Check("popcount of a pattern", Bits.PopCount(0xF0F0F0F0F0F0F0F1UL) == 33);
        }

        private static void CheckWindowsAndSources(Reporter reporter)
        {
            const ulong word = 0xF000000000000001UL;
            reporter.Check("high bits", Windows.Take(word, Window.High, 4) == 0xF);
            reporter.Check("low bits", Windows.Take(word, Window.Low, 4) == 1);
            reporter.Check("high edge bit", Windows.EdgeBit(word, Window.High) == 1);
            reporter.Check("low edge bit", Windows.EdgeBit(word ^ 1, Window.Low) == 0);

            Func<ulong, WordSource> make = seed => new GeneratorWords<SplitMix64>(new SplitMix64(seed));
            var adjacent = new WordStream(new AdjacentSeedWords(make, firstSeed: 10, wordsPerSeed: 2));
            var expectedFirst = new SplitMix64(10);
            var expectedSecond = new SplitMix64(11);
            bool adjacentMatches = adjacent.Next() == expectedFirst.NextUInt64() &&
                adjacent.Next() == expectedFirst.NextUInt64() && adjacent.Next() == expectedSecond.NextUInt64() &&
                adjacent.Next() == expectedSecond.NextUInt64();
            reporter.Check("adjacent seeds give each seed's first words in turn", adjacentMatches);
            bool threw = false;
            try
            {
                new AdjacentSeedWords(make, 1, 3).Fill(new ulong[10], 0, 10);
            }
            catch (ArgumentException)
            {
                threw = true;
            }

            reporter.Check("a fill must hold whole seeds", threw);

            var interleaved = new WordStream(new InterleavedWords(new[] { make(1), make(2) }));
            var one = new SplitMix64(1);
            var two = new SplitMix64(2);
            bool interleavedMatches = interleaved.Next() == one.NextUInt64() &&
                interleaved.Next() == two.NextUInt64() && interleaved.Next() == one.NextUInt64() &&
                interleaved.Next() == two.NextUInt64();
            reporter.Check("interleaved streams take one word each in turn", interleavedMatches);

            var system = new WordStream(new SystemRandomWords(5));
            var again = new WordStream(new SystemRandomWords(5));
            reporter.Check("System.Random words repeat for a seed", system.Next() == again.Next());

            reporter.Check("a verdict under the failure line", new Outcome(1e-7, "").Verdict == "FAIL");
            reporter.Check("a verdict under the suspect line", new Outcome(1e-4, "").Verdict == "SUSPECT");
            reporter.Check("a verdict above both", new Outcome(0.5, "").Verdict == "ok");
            reporter.Check("information has no verdict", Outcome.Information("x").Verdict == "info");

            var even = new List<double>();
            for (int i = 0; i < 200; i++)
            {
                even.Add((i + 0.5) / 200);
            }

            reporter.Check("an even spread is uniform", Calibration.UniformityP(even) > 0.99);
            reporter.Check("a pile at zero is not uniform",
                Calibration.UniformityP(new List<double>(new double[50])) < 1e-6);
        }

        private static int SettledDraws(Outcome outcome)
        {
            string marker = "look fair: ";
            int start = outcome.Detail.IndexOf(marker, StringComparison.Ordinal) + marker.Length;
            int end = outcome.Detail.IndexOf(' ', start);
            return int.Parse(outcome.Detail.Substring(start, end - start));
        }

        private static void CheckControlsAndGenerators(Reporter reporter, TextWriter output)
        {
            var failsEach = new Dictionary<string, string[]>
            {
                ["RANDU"] = new[] { "serial, high bits", "birthday spacings 2^32, high bits",
                    "binary rank 32x32, high bits" },
                ["MINSTD"] = new[] { "birthday spacings 2^28, high bits", "birthday spacings 2^32, high bits",
                    "binary rank 32x32, high bits" },
                ["LCG64"] = new[] { "serial, low bits", "linear complexity, low bits" },
                ["xorshift64"] = new[] { "linear complexity, high bits", "linear complexity, low bits",
                    "binary rank 64x64" },
                ["xoshiro256+"] = new[] { "linear complexity, low bits" },
                ["System.Random"] = new[] { "first 8 draws of each seed, high bits",
                    "first draw of consecutive seeds, high bits" },
            };
            List<Subject> subjects = Subjects.All();
            output.WriteLine("Running the quick battery on " + subjects.Count + " subjects.");
            List<NamedOutcome>[] outcomes = Battery.RunAll(subjects, new Scale(Program.QuickDivisor),
                Program.DefaultSeed, subject => { });
            for (int i = 0; i < subjects.Count; i++)
            {
                Subject subject = subjects[i];
                if (subject.IsControl)
                {
                    foreach (string name in failsEach[subject.Name])
                    {
                        NamedOutcome found = outcomes[i].FirstOrDefault(outcome => outcome.Name == name);
                        reporter.Check(subject.Name + " fails " + name,
                            found != null && found.Outcome.Verdict == "FAIL",
                            found == null ? "no such check" : found.Outcome.PText);
                    }

                    continue;
                }

                foreach (NamedOutcome outcome in outcomes[i])
                {
                    if (outcome.ExpectsFailure)
                    {
                        reporter.Check(subject.Name + " fails its control '" + outcome.Name + "'",
                            outcome.Outcome.Verdict == "FAIL");
                    }
                    else
                    {
                        reporter.Check(subject.Name + " passes " + outcome.Name, outcome.Outcome.Verdict != "FAIL",
                            outcome.Outcome.PText);
                    }
                }
            }

            int splitMixIndex = subjects.FindIndex(subject => subject.Name == "SplitMix64");
            int xoshiroIndex = subjects.FindIndex(subject => subject.Name == "xoshiro256++");
            Func<int, int> settled = index => SettledDraws(outcomes[index].First(
                outcome => outcome.Name == "recovery from single-bit states").Outcome);
            reporter.Check("SplitMix64 mixes a single-bit state at once", settled(splitMixIndex) == 0);
            reporter.Check("xoshiro256++ needs some draws to mix a single-bit state", settled(xoshiroIndex) > 0);
        }

        private static void CheckRanges(Reporter reporter)
        {
            reporter.Check("gap below 1.5", RangeChecks.GapBelow(1.5, RangeChecks.DoubleMantissaBits) ==
                Math.Pow(2, -52));
            reporter.Check("gap below a power of two is the finer one", RangeChecks.GapBelow(1,
                RangeChecks.DoubleMantissaBits) == Math.Pow(2, -53));
            reporter.Check("gap below 7.25", RangeChecks.GapBelow(7.25, RangeChecks.DoubleMantissaBits) ==
                Math.Pow(2, -50));
            reporter.Check("gap below a negative end", RangeChecks.GapBelow(-1000, RangeChecks.DoubleMantissaBits) ==
                Math.Pow(2, -43));
            reporter.Check("float gap below 1.5", RangeChecks.GapBelow(1.5, RangeChecks.FloatMantissaBits) ==
                Math.Pow(2, -23));
            reporter.Check("float gap below 2", RangeChecks.GapBelow(2, RangeChecks.FloatMantissaBits) ==
                Math.Pow(2, -23));

            var largest = new LargestWords();
            reporter.Check("the largest double is one step under 1", largest.NextDouble() == 1 - Math.Pow(2, -53));
            reporter.Check("the largest float is one step under 1", largest.NextFloat() == 1 - (float)Math.Pow(2, -24));

            var generator = new Sfc64(1);
            double gap = Math.Pow(2, -52);
            NamedOutcome even = RangeChecks.OnGrid("even", 1, 2, mantissaBits: RangeChecks.DoubleMantissaBits,
                draws: 1 << 20, draw: () => 1 + generator.NextUInt64(1UL << 52) * gap);
            reporter.Check("draws uniform over the cells pass", even.Outcome.Verdict == "ok", even.Outcome.PText);
            NamedOutcome negative = RangeChecks.OnGrid("negative", -1.5, -1,
                mantissaBits: RangeChecks.DoubleMantissaBits, draws: 1 << 20,
                draw: () => -1.5 + generator.NextUInt64(1UL << 51) * gap);
            reporter.Check("draws uniform over negative cells pass", negative.Outcome.Verdict == "ok",
                negative.Outcome.PText);
            NamedOutcome evenCells = RangeChecks.OnGrid("even cells", 1, 2,
                mantissaBits: RangeChecks.DoubleMantissaBits, draws: 1 << 20,
                draw: () => 1 + generator.NextUInt64(1UL << 51) * 2 * gap);
            reporter.Check("draws that never land on an odd cell fail", evenCells.Outcome.Verdict == "FAIL");
        }

        public static int Run(TextWriter output)
        {
            var reporter = new Reporter(output);
            CheckStatistics(reporter);
            CheckRanges(reporter);
            CheckLinearAlgebra(reporter);
            CheckWindowsAndSources(reporter);
            CheckControlsAndGenerators(reporter, output);
            output.WriteLine(reporter.Failures == 0
                ? "All " + reporter.Checks + " self-checks passed."
                : reporter.Failures + " of " + reporter.Checks + " self-checks FAILED.");
            return reporter.Failures == 0 ? 0 : 1;
        }
    }
}