Orthogonal fractional factorial experiment design

Program.cs

using System.Diagnostics;
using System.Runtime.CompilerServices;
using System.Runtime.InteropServices;

namespace ExperimentDesign
{
    internal class Program
    {
        static void Main(string[] args)
        {
            Console.WriteLine("DOE orthogonal plan generator.");
            Console.WriteLine("Input number of levels for each factors (one per line):");

            List<uint> levels = [];
            while (true)
            {
                var line = Console.ReadLine();
                if (!uint.TryParse(line, out var level) || level == 0)
                {
                    if (!string.IsNullOrWhiteSpace(line))
                    {
                        Console.WriteLine("Invalid number, assuming end of input");
                    }
                    break;
                }
                levels.Add(level);
            }

            if (levels.Count == 0)
            {
                Console.WriteLine("No input provided.");
                return;
            }

            var plan = new Plan(levels);
            plan.Generate();

            var result = new uint[levels.Count];
            for (uint i = 0; i < plan.SampleCount; ++i)
            {
                plan.GetSamples(i, result);
                foreach (var v in result)
                {
                    Console.Write(v);
                    Console.Write('\t');
                }
                Console.WriteLine();
            }
        }

        internal abstract class OrthogonalArray
        {
            protected readonly uint _base, _power, _size;

            public abstract void SetColumns(uint n);
            public abstract uint GetCell(uint row, uint column);

            public static OrthogonalArray Create(uint nbase, uint power, uint size)
            {
                Debug.Assert(Math.Pow(nbase, power) == size);
                ArgumentOutOfRangeException.ThrowIfGreaterThan(nbase, 255u);
                if (power == 1)
                {
                    return new SingleColumnArray(nbase, power, size);
                }
                else
                {
                    return new GeneralArray(nbase, power, size);
                }
            }

            private OrthogonalArray(uint nbase, uint power, uint size)
            {
                _base = nbase;
                _power = power;
                _size = size;
            }

            private sealed class SingleColumnArray : OrthogonalArray
            {
                public SingleColumnArray(uint nbase, uint power, uint size)
                    : base(nbase, power, size)
                {
                    Debug.Assert(power == 1);
                }

                public override void SetColumns(uint n)
                {
                    ArgumentOutOfRangeException.ThrowIfNotEqual(n, 1u);
                }

                public override uint GetCell(uint row, uint column)
                {
                    ArgumentOutOfRangeException.ThrowIfNotEqual(column, 0u);
                    ArgumentOutOfRangeException.ThrowIfGreaterThanOrEqual(row, _base);
                    return row;
                }
            }

            /*
             * Example for 3^3=27 trials experiments, which can holds 13 factors, as 13*(3-1)=27-1.
             * Three elementary factors are called A, B, and C.
             * Other factors are combinations of them.
             * 
             * C
             * B, B+C, B+2C
             * A, A+B, A+2B, A+C, A+B+C, A+2B+C, A+2C, A+B+2C, A+2B+2C
             * 
             * One can confirm any pair among these factors are orthogonal.
             * The columns are reordered vertically: C, B, A, B+C, A+B, B+2C, A+2B, A+C, ...
             * This can avoid adjacent factors to be correlated.
             * 
             * In the implementation below, A+B+2C is represented as [2,1,1].
             * Representations for all columns are concatenated and saved in the list for later use.
             */
            private sealed class GeneralArray : OrthogonalArray
            {
                private readonly List<byte> _columns = [];
                private readonly uint _maxColumns;
                private uint _columnCount = 0;

                public GeneralArray(uint nbase, uint power, uint size)
                    : base(nbase, power, size)
                {
                    _maxColumns = (size - 1) / (nbase - 1);
                }

                public override void SetColumns(uint n)
                {
                    ArgumentOutOfRangeException.ThrowIfGreaterThan(n, _maxColumns);
                    _ = checked((int)(n * _power));
                    _columns.Clear();

                    //Implement the reordering of columns.
                    var power = _power;
                    var nbase = _base;
                    var coeff = new byte[power + 1];
                    var startingLine = 0;
                    while ((uint)_columns.Count < n * power)
                    {
                        for (int line = startingLine; line < power; ++line)
                        {
                            for (int i = startingLine; i < line; ++i)
                            {
                                _columns.Add(0);
                            }
                            for (int i = 0; i < startingLine; ++i)
                            {
                                _columns.Add(coeff[i]);
                            }
                            _columns.Add(1);
                            for (int i = line + 1; i < power; ++i)
                            {
                                _columns.Add(0);
                            }
                        }

                        int inc = 0;
                        while (true)
                        {
                            if (++coeff[inc] < nbase)
                            {
                                break;
                            }
                            coeff[inc] = 0;
                            inc += 1;
                        }
                        if (inc >= startingLine)
                        {
                            startingLine = inc + 1;
                        }
                    }
                    _columnCount = n;
                }

                public override uint GetCell(uint row, uint column)
                {
                    ArgumentOutOfRangeException.ThrowIfGreaterThanOrEqual(row, _size);
                    ArgumentOutOfRangeException.ThrowIfGreaterThanOrEqual(column, _columnCount);

                    var coeff = CollectionsMarshal.AsSpan(_columns)[(int)(column * _power)..];
                    uint sum = 0;
                    for (uint i = 0; i < _power; ++i)
                    {
                        (row, var r) = Math.DivRem(row, _base);
                        sum = (sum + r * coeff[(int)i]) % _base;
                    }

                    return sum;
                }
            }
        }

        private sealed class PrimeSeqence
        {
            private static readonly uint[] _values =
            [
                2, 3, 5, 7, 11, 13, 17, 19,
                23, 29, 31, 37, 41, 43, 47,
            ];

            public uint Get(int index)
            {
                return _values[index];
            }
        }

        private sealed class Plan
        {
            private sealed class PseudoFactorGroup
            {
                private uint _count = 0;
                private uint _samples = 0;
                private OrthogonalArray? _oa;

                public uint Level { get; }
                public uint SampleCount => _samples == 0 ? throw new InvalidOperationException() : _samples;

                public PseudoFactorGroup(uint level)
                {
                    Debug.Assert(level >= 2);
                    Level = level;
                }

                public uint AddFactor()
                {
                    Debug.Assert(_samples == 0);
                    return _count++;
                }

                public void Generate()
                {
                    Debug.Assert(_samples == 0);
                    var level = Level;
                    var dof = (level - 1) * _count;
                    uint power = 0;
                    uint exp = 1;
                    while (exp <= dof)
                    {
                        power += 1;
                        exp = checked(exp * level);
                    }
                    _samples = exp;
                    _oa = OrthogonalArray.Create(level, power, exp);
                    _oa.SetColumns(_count);
                }

                public uint Get(uint index, uint column)
                {
                    var oa = _oa ?? throw new InvalidOperationException();
                    return oa.GetCell(index, column);
                }
            }

            private readonly struct FactorMapping
            {
                public int PseudoFactorIndex { get; init; }
                public uint ColumnIndex { get; init; }
            }

            private readonly PrimeSeqence _primes = new();

            private readonly List<uint> _levels;
            private readonly List<PseudoFactorGroup> _pseudoFactors = [];
            private readonly Dictionary<uint, int> _pseudoFactorsByLevel = [];
            private readonly List<List<FactorMapping>> _factorMappings = [];

            private uint[] _pseudoFactorIndicesBuffer = null!;
            private uint _samples = 0;

            public uint SampleCount => _samples == 0 ? throw new InvalidOperationException() : _samples;

            public Plan(List<uint> levels)
            {
                _levels = levels;
            }

            public void Generate()
            {
                SplitFactors();
                PreparePseudoFactors();
            }

            private void SplitFactors()
            {
                foreach (uint l in _levels)
                {
                    uint level = l;
                    var mapping = new List<FactorMapping>();
                    _factorMappings.Add(mapping);

                    if (level > 1)
                    {
                        FindAndAddPseudoFactor(ref level, 2, mapping);
                    }
                    if (level > 1)
                    {
                        FindAndAddPseudoFactor(ref level, 3, mapping);
                    }
                    if (level > 1)
                    {
                        FindAndAddPseudoFactor(ref level, 5, mapping);
                    }
                    if (level > 1)
                    {
                        FindAndAddPseudoFactor(ref level, 7, mapping);
                    }
                    int pindex = 4;
                    while (level > 1)
                    {
                        FindAndAddPseudoFactor(ref level, _primes.Get(pindex++), mapping);
                    }
                }
            }

            [MethodImpl(MethodImplOptions.AggressiveInlining)]
            private void FindAndAddPseudoFactor(ref uint level, uint pf, List<FactorMapping> mapping)
            {
                while (level > 1)
                {
                    var (q, r) = Math.DivRem(level, pf);
                    if (r != 0)
                    {
                        return;
                    }

                    if (!_pseudoFactorsByLevel.TryGetValue(pf, out var index))
                    {
                        index = _pseudoFactors.Count;
                        _pseudoFactorsByLevel[pf] = index;
                        _pseudoFactors.Add(new PseudoFactorGroup(pf));
                    }
                    var pfo = _pseudoFactors[index];
                    var col = pfo.AddFactor();
                    mapping.Add(new()
                    {
                        PseudoFactorIndex = index,
                        ColumnIndex = col,
                    });
                    level = q;
                }
            }

            private void PreparePseudoFactors()
            {
                uint samples = 1;
                foreach (var pf in _pseudoFactors)
                {
                    pf.Generate();
                    samples = checked(samples * pf.SampleCount);
                }
                _samples = samples;

                _pseudoFactorIndicesBuffer = new uint[_pseudoFactors.Count];
            }

            public void GetSamples(uint index, Span<uint> result)
            {
                if (result.Length != _levels.Count)
                {
                    throw new ArgumentException("length of result Span must be equal to factor count");
                }

                for (int i = 0; i < _pseudoFactors.Count; ++i)
                {
                    var pf = _pseudoFactors[i];
                    (index, _pseudoFactorIndicesBuffer[i]) = Math.DivRem(index, pf.SampleCount);
                }

                for (int i = 0; i < _levels.Count; ++i)
                {
                    var mapping = _factorMappings[i];
                    uint multiplier = 1;
                    uint sum = 0;
                    foreach (var m in mapping)
                    {
                        var pf = _pseudoFactors[m.PseudoFactorIndex];
                        var pfv = pf.Get(_pseudoFactorIndicesBuffer[m.PseudoFactorIndex], m.ColumnIndex);
                        sum += pfv * multiplier;
                        multiplier *= pf.Level;
                    }
                    result[i] = sum;
                }
            }
        }
    }
}