Skip to content

Instantly share code, notes, and snippets.

@acaly
Created August 3, 2024 15:23
Show Gist options
  • Save acaly/127af291d0f5747fea1cc29aa4457cac to your computer and use it in GitHub Desktop.
Save acaly/127af291d0f5747fea1cc29aa4457cac to your computer and use it in GitHub Desktop.
Orthogonal fractional factorial experiment design
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;
}
}
}
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment